home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / phigs / ptk.lha / ptk / source / library / tran.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-01  |  85.0 KB  |  2,893 lines
  1. /*----------------------------------------------------------------------------
  2.  
  3.  Module name: transformations utility package. 
  4.  
  5.  Authors: W.T. Hewitt, K.M. Singleton.
  6.  
  7.  Function: contains functions for performing two-dimensional and 
  8.  three-dimensional transformations. 
  9.  
  10.  Dependencies: machine.h, transformtypes.h.
  11.  
  12.  Internal function list: validbounds.
  13.  
  14.  External function list: ptk_equal, ptk_point, ptk_point3, ptk_vector,
  15.  ptk_vector3, ptk_vec3topt3, ptk_pt3tovec3, ptk_limit, ptk_limit3,
  16.  ptk_dotv, ptk_crossv,
  17.  ptk_nullv, ptk_modv, ptk_unitv, ptk_scalev, ptk_subv, ptk_addv,
  18.  ptk_unitmatrix3, ptk_transposematrix3, ptk_multiplymatrix3, 
  19.  ptk_concatenatematrix3,
  20.  ptk_shift3, ptk_scale3, ptk_rotatecs3, ptk_rotate3, ptk_shear3,
  21.  ptk_rotatevv3, ptk_rotateline3, ptk_pt3topt4, ptk_pt4topt3, ptk_transform4,
  22.  ptk_transform3, ptk_matrixtomatrix3, ptk_outputmatrix3, ptk_box3tobox3, 
  23.  ptk_accumulatetran3, ptk_evalvieworientation3,
  24.  ptk_evalviewmapping3, ptk_stackmatrix3, ptk_unstackmatrix3, 
  25.  ptk_examinestackmatrix3,
  26.  ptk_3ptto3pt, ptk_0to3pt, ptk_oto3pt, ptk_invertmatrix3.
  27.  
  28.  Hashtables used: none.
  29.  
  30.  Modification history: (Version), (Date), (Name), (Description).
  31.  
  32.  1.0,  2nd February 1986,  W.T. Hewitt,  First version.
  33.  
  34.  1.1,  1st April 1986,  K.M. Singleton, W.T. Hewitt, Added viewing 
  35.  and other utility routines.
  36.  
  37.  1.2,  8th June 1987,  K.M. Wyrwas,  Added routine to invert a matrix.
  38.  
  39.  1.3,  8th July 1988,  S.Larkin,  Modified to work with VAX PHIGS$.
  40.  
  41.  1.4,  30th October 1990,  W.T. Hewitt,   Fixed to work with VAX PHIGS$.
  42.  
  43.  2.0,  4th January 1991,  G. Williams,  Converted from Pascal to C.
  44.  
  45.  2.1,  13th February 1991, G. Williams, Added functions - ptk_vector3,
  46.  ptk_vec3topt3 and ptk_pt3tovec3.
  47.  
  48.  2.2, 3rd May 1991, G. Williams, Added function - ptk_3x3to4x4.
  49.  
  50.  3.0, June 1992, G. Williams, Converted to ISO PHIGS C.
  51.  
  52. ----------------------------------------------------------------------------*/
  53.  
  54. #include <stdio.h>
  55. #include <math.h>
  56. #ifdef SUN
  57. #include <malloc.h>
  58. #include <memory.h>
  59. #endif
  60. #include "machine.h"
  61. #include "ptktype.h"
  62. #include "trantype.h"
  63.  
  64. /* Pointer records */
  65.  
  66. typedef struct ptksstack3 
  67. {
  68.   Pmatrix3 ptkamat;
  69.   struct ptksstack3 *ptkpnext;
  70. } ptksstack3;
  71.  
  72. typedef struct ptksstack
  73. {
  74.   Pmatrix ptkamat;
  75.   struct ptksstack *ptkpnext;
  76. } ptksstack;
  77.  
  78. #define ptkcdtor         (3.14159263 / 180.0)
  79.  
  80. static ptksstack3 *listhead3 = NULL, *newone3 = NULL;
  81. static ptksstack *listhead = NULL, *newone = NULL;
  82.  
  83. /* global variable for ptkcpceps */
  84.  
  85. Pfloat ptkveps = ptkcpceps;
  86.  
  87. /*--------------------------------------------------------------------------*/
  88.  
  89. /*function:external*/
  90. extern ptkboolean ptk_equal(C(Pfloat) one, C(Pfloat) two)
  91. PreANSI(Pfloat one)
  92. PreANSI(Pfloat two)
  93. /*
  94. ** \parambegin
  95. ** \param{Pfloat}{one}{floating point number}{IN}
  96. ** \param{Pfloat}{two}{floating point number}{IN}
  97. ** \paramend
  98. ** \blurb{This function returns TRUE if \pardesc{one} and \pardesc{two}
  99. ** are equal, or their difference is less than the global 
  100. ** constant tolerance \pardesc{ptkveps}. This is a global 
  101. ** variable of type \pardesc{Pfloat} which may be changed by the application.
  102. ** The default value of \pardesc{ptkveps} is 1.0e-7.}
  103. */
  104. {
  105.   Pfloat x, y;
  106.  
  107.   x = (Pfloat)one;
  108.   y = (Pfloat)two;
  109.   return (fabs(x - y) <= ptkveps);
  110. }  /* ptk_equal */
  111.  
  112. /*--------------------------------------------------------------------------*/
  113.  
  114. /*function:external*/
  115. extern Ppoint ptk_point(C(Pfloat) x, C(Pfloat) y)
  116. PreANSI(Pfloat x)
  117. PreANSI(Pfloat y)
  118. /*
  119. ** \parambegin
  120. ** \param{Pfloat}{x}{x coordinate}{IN}
  121. ** \param{Pfloat}{y}{y coordinate}{IN}
  122. ** \paramend
  123. ** \blurb{This function returns a structure representing the 
  124. ** 2D point \pardesc{(x,y)}.}
  125. */
  126. {
  127.   Ppoint temp;
  128.  
  129.   temp.x = x;
  130.   temp.y = y;
  131.   return temp;
  132. }  /* ptk_point */
  133.  
  134. /*--------------------------------------------------------------------------*/
  135.  
  136. /*function:external*/
  137. extern Ppoint3 ptk_point3(C(Pfloat) x, C(Pfloat) y, C(Pfloat) z)
  138. PreANSI(Pfloat x)
  139. PreANSI(Pfloat y)
  140. PreANSI(Pfloat z)
  141. /*
  142. ** \parambegin
  143. ** \param{Pfloat}{x}{x coordinate}{IN}
  144. ** \param{Pfloat}{y}{y coordinate}{IN}
  145. ** \param{Pfloat}{z}{z coordinate}{IN}
  146. ** \paramend
  147. ** \blurb{This function returns a structure representing
  148. ** the 3D point \pardesc{(x,y,z)}.}
  149. */
  150. {
  151.   Ppoint3 temp;
  152.  
  153.   temp.x = x;
  154.   temp.y = y;
  155.   temp.z = z;
  156.   return temp;
  157. }  /* ptk_point3 */
  158.  
  159. /*--------------------------------------------------------------------------*/
  160.  
  161. /*function:external*/
  162. extern Pvec ptk_vector(C(Pfloat) x, C(Pfloat) y)
  163. PreANSI(Pfloat x)
  164. PreANSI(Pfloat y)
  165. /*
  166. ** \parambegin
  167. ** \param{Pfloat}{x}{x coordinate}{IN}
  168. ** \param{Pfloat}{y}{y coordinate}{IN}
  169. ** \paramend
  170. ** \blurb{This function returns a structure representing
  171. ** the 2D vector \pardesc{(x,y)}.}
  172. */
  173. {
  174.   Pvec temp;
  175.  
  176.   temp.delta_x = x;
  177.   temp.delta_y = y;
  178.   return temp;
  179. }  /* ptk_vector */
  180.  
  181. /*--------------------------------------------------------------------------*/
  182.  
  183. /*function:external*/
  184. extern Pvec3 ptk_vector3(C(Pfloat) x, C(Pfloat) y, C(Pfloat) z)
  185. PreANSI(Pfloat x)
  186. PreANSI(Pfloat y)
  187. PreANSI(Pfloat z)
  188. /*
  189. ** \parambegin
  190. ** \param{Pfloat}{x}{x coordinate}{IN}
  191. ** \param{Pfloat}{y}{y coordinate}{IN}
  192. ** \param{Pfloat}{z}{z coordinate}{IN}
  193. ** \paramend
  194. ** \blurb{This function returns a structure representing
  195. ** the 3D vector \pardesc{(x,y,z)}.}
  196. */
  197. {
  198.   Pvec3 temp;
  199.  
  200.   temp.delta_x = x;
  201.   temp.delta_y = y;
  202.   temp.delta_z = z;
  203.   return temp;
  204. }  /* ptk_vector3 */
  205.  
  206. /*--------------------------------------------------------------------------*/
  207.  
  208. /*function:external*/
  209. extern Plimit ptk_limit(C(Pfloat) xmin, C(Pfloat) xmax, C(Pfloat) ymin,
  210.                         C(Pfloat) ymax)
  211. PreANSI(Pfloat xmin)
  212. PreANSI(Pfloat xmax)
  213. PreANSI(Pfloat ymin)
  214. PreANSI(Pfloat ymax)
  215. /*
  216. ** \parambegin
  217. ** \param{Pfloat}{xmin}{minimum x coordinate}{IN}
  218. ** \param{Pfloat}{xmax}{maximum x coordinate}{IN}
  219. ** \param{Pfloat}{ymin}{minimum y coordinate}{IN}
  220. ** \param{Pfloat}{ymax}{maximum y coordinate}{IN}
  221. ** \paramend
  222. ** \blurb{This function returns a structure representing the 2D
  223. ** rectangle whose corners are defined by \pardesc{xmin, xmax, ymin, ymax}.}
  224. */
  225. {
  226.   Plimit temp;
  227.  
  228.   temp.x_min = xmin;
  229.   temp.x_max = xmax;
  230.   temp.y_min = ymin;
  231.   temp.y_max = ymax;
  232.   return temp;
  233. }  /* ptk_limit */
  234.  
  235. /*--------------------------------------------------------------------------*/
  236.  
  237. /*function:external*/
  238. extern Plimit3 ptk_limit3(C(Pfloat) xmin, C(Pfloat) xmax, C(Pfloat) ymin,
  239.                           C(Pfloat) ymax, C(Pfloat) zmin, C(Pfloat) zmax)
  240. PreANSI(Pfloat xmin)
  241. PreANSI(Pfloat xmax)
  242. PreANSI(Pfloat ymin)
  243. PreANSI(Pfloat ymax)
  244. PreANSI(Pfloat zmin)
  245. PreANSI(Pfloat zmax)
  246. /*
  247. ** \parambegin
  248. ** \param{Pfloat}{xmin}{minimum x coordinate}{IN}
  249. ** \param{Pfloat}{xmax}{maximum x coordinate}{IN}
  250. ** \param{Pfloat}{ymin}{minimum y coordinate}{IN}
  251. ** \param{Pfloat}{ymax}{maximum y coordinate}{IN}
  252. ** \param{Pfloat}{zmin}{minimum z coordinate}{IN}
  253. ** \param{Pfloat}{zmax}{maximum z coordinate}{IN}
  254. ** \paramend
  255. ** \blurb{This function returns a structure representing the 
  256. ** 3D volume defined by \pardesc{xmin, xmax, ymin, ymax, zmin, zmax}.}
  257. */
  258. {
  259.   Plimit3 temp;
  260.  
  261.   temp.x_min = xmin;
  262.   temp.x_max = xmax;
  263.   temp.y_min = ymin;
  264.   temp.y_max = ymax;
  265.   temp.z_min = zmin;
  266.   temp.z_max = zmax;
  267.   return temp;
  268. }  /* ptk_limit3 */
  269.  
  270. /*--------------------------------------------------------------------------*/
  271.  
  272. /*function:external*/
  273. extern Ppoint3 ptk_vec3topt3(C(Pvec3 *) vec)
  274. PreANSI(Pvec3 *vec)
  275. /*
  276. ** \parambegin
  277. ** \param{Pvec3 *}{vec}{3D vector}{IN}
  278. ** \paramend
  279. ** \blurb{This function convert the 3D vector \pardesc{vec} to a 3D point.}
  280. */
  281. {
  282.   Ppoint3 temp;
  283.  
  284.   temp = ptk_point3(vec->delta_x, vec->delta_y, vec->delta_z);
  285.   return temp;
  286. }  /* ptk_vec3topt3 */
  287.  
  288. /*--------------------------------------------------------------------------*/
  289.  
  290. /*function:external*/
  291. extern Pvec3 ptk_pt3tovec3(C(Ppoint3 *) pt)
  292. PreANSI(Ppoint3 *pt)
  293. /*
  294. ** \parambegin
  295. ** \param{Ppoint3 *}{pt}{3D point}{IN}
  296. ** \paramend
  297. ** \blurb{This function converts the 3D point \pardesc{pt} to a 3D vector.}
  298. */
  299. {
  300.   Pvec3 temp;
  301.  
  302.   temp = ptk_vector3(pt->x, pt->y, pt->z);
  303.   return temp;
  304. }  /* ptk_pt3tovec3 */
  305.  
  306. /*--------------------------------------------------------------------------*/
  307.  
  308. /*function:external*/
  309. extern Pfloat ptk_dotv3(C(Ppoint3 *) v1, C(Ppoint3 *) v2)
  310. PreANSI(Ppoint3 *v1)
  311. PreANSI(Ppoint3 *v2)
  312. /*
  313. ** \parambegin
  314. ** \param{Ppoint3 *}{v1}{3D vector}{IN}
  315. ** \param{Ppoint3 *}{v2}{3D vector}{IN}
  316. ** \paramend
  317. ** \blurb{This function evaluates the dot product of the
  318. **  two 3D vectors \pardesc{v1} and 
  319. ** \pardesc{v2}, returning it as the value of the function.}
  320. */
  321. {
  322.   return (v1->x * v2->x + v1->y * v2->y + v1->z * v2->z);
  323. }  /* ptk_dotv3 */
  324.  
  325. /*--------------------------------------------------------------------------*/
  326.  
  327. /*function:external*/
  328. extern Pfloat ptk_dotv(C(Ppoint *) v1, C(Ppoint *) v2)
  329. PreANSI(Ppoint *v1)
  330. PreANSI(Ppoint *v2)
  331. /*
  332. ** \parambegin
  333. ** \param{Ppoint *}{v1}{2D vector}{IN}
  334. ** \param{Ppoint *}{v2}{2D vector}{IN}
  335. ** \paramend
  336. ** \blurb{Evaluates the dot product of the two 2D vectors \pardesc{v1} and
  337. ** \pardesc{v2}, returning it as the value of the function.}
  338. */
  339. {
  340.   return (v1->x * v2->x + v1->y * v2->y);
  341. }  /* ptk_dotv */
  342.  
  343. /*--------------------------------------------------------------------------*/
  344.  
  345. /*function:external*/
  346. extern Ppoint3 ptk_crossv3(C(Ppoint3 *) v1, C(Ppoint3 *) v2)
  347. PreANSI(Ppoint3 *v1)
  348. PreANSI(Ppoint3 *v2)
  349. /*
  350. ** \parambegin
  351. ** \param{Ppoint3 *}{v1}{3D vector}{IN}
  352. ** \param{Ppoint3 *}{v2}{3D vector}{IN}
  353. ** \paramend
  354. ** \blurb{This function evaluates the cross product of
  355. ** the two 3D vectors \pardesc{v1} and
  356. ** \pardesc{v2}, returning the new vector as
  357. ** the function result. Since a local copy is made statements such as
  358. ** {\tt v2 == ptk\_crossv(v1, v2)}
  359. ** will produce the correct answer.}
  360. */
  361. {
  362.   Ppoint3 temp;
  363.  
  364.   temp.x = v1->y * v2->z - v1->z * v2->y;
  365.   temp.y = v1->z * v2->x - v1->x * v2->z;
  366.   temp.z = v1->x * v2->y - v1->y * v2->x;
  367.   return temp;
  368. }  /* ptk_crossv3 */
  369.  
  370. /*--------------------------------------------------------------------------*/
  371.  
  372. /*function:external*/
  373. extern ptkboolean ptk_nullv3(C(Ppoint3 *) vec)
  374. PreANSI(Ppoint3 *vec)
  375. /*
  376. ** \parambegin
  377. ** \param{Ppoint3 *}{vec}{3D vector}{IN}
  378. ** \paramend
  379. ** \blurb{This function returns \pardesc{TRUE} if the modulus
  380. ** of the 3D vector\pardesc{vec} 
  381. ** is less than the global tolerance \pardesc{ptkpceps}, otherwise 
  382. ** \pardesc{FALSE}.}
  383. */
  384. {
  385.   return (ptk_equal(vec->x, 0.0) && ptk_equal(vec->y, 0.0) &&
  386.       ptk_equal(vec->z, 0.0));
  387. }  /* ptk_nullv3 */
  388.  
  389. /*--------------------------------------------------------------------------*/
  390.  
  391. /*function:external*/
  392. extern ptkboolean ptk_nullv(C(Ppoint *) vec)
  393. PreANSI(Ppoint *vec)
  394. /*
  395. ** \parambegin
  396. ** \param{Ppoint *}{vec}{2D vector}{IN}
  397. ** \paramend
  398. ** \blurb{This function returns \pardesc{TRUE} if the
  399. ** modulus of the 2D vector \pardesc{vec} 
  400. ** is less than the global tolerance \pardesc{ptkpceps}, otherwise 
  401. ** \pardesc{FALSE}.}
  402. */
  403. {
  404.   return (ptk_equal(vec->x, 0.0) && ptk_equal(vec->y, 0.0));
  405. }  /* ptk_nullv */
  406.  
  407. /*--------------------------------------------------------------------------*/
  408.  
  409. /*function:external*/
  410. extern Pfloat ptk_modv3(C(Ppoint3 *) vec)
  411. PreANSI(Ppoint3 *vec)
  412. /*
  413. ** \parambegin
  414. ** \param{Ppoint3 *}{vec}{3D vector}{IN}
  415. ** \paramend
  416. ** \blurb{Returns the modulus of the vector \pardesc{vec}.}
  417. */
  418. {
  419.   return sqrt(vec->x * vec->x + vec->y * vec->y + vec->z * vec->z);
  420. }  /* ptk_modv3 */
  421.  
  422. /*-------------------------------------------------------------------------*/
  423.  
  424. /*function:external*/
  425. extern Pfloat ptk_modv(C(Ppoint *) vec)
  426. PreANSI(Ppoint *vec)
  427. /*
  428. ** \parambegin
  429. ** \param{Ppoint *}{vec}{2D vector}{IN}
  430. ** \paramend
  431. ** \blurb{This function returns the modulus of the 2D vector \pardesc{vec}.}
  432. */
  433. {
  434.   return sqrt(vec->x * vec->x + vec->y * vec->y);
  435. }  /* ptk_modv */
  436.  
  437. /*-------------------------------------------------------------------------*/
  438.  
  439. /*function:external*/
  440. extern Ppoint3 ptk_unitv3(C(Ppoint3 *) vec)
  441. PreANSI(Ppoint3 *vec)
  442. /*
  443. ** \parambegin
  444. ** \param{Ppoint3 *}{vec}{3D vector}{IN}
  445. ** \paramend
  446. ** \blurb{This function generates and returns a unit
  447. ** vector from the supplied 3D vector 
  448. ** \pardesc{vec}.}
  449. */
  450. {
  451.   Pfloat modu;
  452.   Ppoint3 temp;
  453.  
  454.   modu = ptk_modv3(vec);
  455.   if (!ptk_equal(modu, 0.0)) 
  456.   {
  457.     temp = ptk_point3(vec->x / modu, vec->y / modu, vec->z / modu);
  458.     return temp;
  459.   } 
  460.   else 
  461.   {
  462.     temp = ptk_point3(0.0, 0.0, 0.0);
  463.     return temp;
  464.   }
  465. }  /* ptk_unitv3 */
  466.  
  467. /*--------------------------------------------------------------------------*/
  468.  
  469. /*function:external*/
  470. extern Ppoint ptk_unitv(C(Ppoint *) vec)
  471. PreANSI(Ppoint *vec)
  472. /*
  473. ** \parambegin
  474. ** \param{Ppoint *}{vec}{2D vector}{IN}
  475. ** \paramend
  476. ** \blurb{This function generates and returns a unit vector from
  477. ** the supplied 2D vector 
  478. ** \pardesc{vec}.}
  479. */
  480. {
  481.   Pfloat modu;
  482.   Ppoint temp;
  483.  
  484.   modu = ptk_modv(vec);
  485.   if (!ptk_equal(modu, 0.0)) 
  486.   {
  487.     temp = ptk_point(vec->x / modu, vec->y / modu);
  488.     return temp;
  489.   } 
  490.   else 
  491.   {
  492.     temp = ptk_point(0.0, 0.0);
  493.     return temp;
  494.   }
  495. }  /* ptk_unitv */
  496.  
  497. /*--------------------------------------------------------------------------*/
  498.  
  499. /*function:external*/
  500. extern Ppoint3 ptk_scalev3(C(Ppoint3 *) vec, C(Pfloat) scale)
  501. PreANSI(Ppoint3 *vec)
  502. PreANSI(Pfloat scale)
  503. /*
  504. ** \parambegin
  505. ** \param{Ppoint3 *}{vec}{3D vector}{IN}
  506. ** \param{Pfloat}{scale}{scale factor}{IN}
  507. ** \paramend
  508. ** \blurb{This function multiplies the 3D vector \pardesc{v}
  509. **  by the  scalar \pardesc{s} and
  510. ** returns the result.}
  511. */
  512. {
  513.   Ppoint3 temp;
  514.  
  515.   temp = ptk_point3(vec->x * scale, vec->y * scale, vec->z * scale);
  516.   return temp;
  517. }  /* ptk_scalev3 */
  518.  
  519. /*--------------------------------------------------------------------------*/
  520.  
  521. /*function:external*/
  522. extern Ppoint ptk_scalev(C(Ppoint *) vec, C(Pfloat) scale)
  523. PreANSI(Ppoint *vec)
  524. PreANSI(Pfloat scale)
  525. /*
  526. ** \parambegin
  527. ** \param{Ppoint *}{vec}{2D vector}{IN}
  528. ** \param{Pfloat}{scale}{scale factor}{IN}
  529. ** \paramend
  530. ** \blurb{This function multiplies the 2D vector
  531. ** \pardesc{v} by the scalar \pardesc{s} and
  532. ** returns the result.}
  533. */
  534. {
  535.   Ppoint temp;
  536.  
  537.   temp = ptk_point(vec->x * scale, vec->y * scale);
  538.   return temp;
  539. }  /* ptk_scalev */
  540.  
  541. /*--------------------------------------------------------------------------*/
  542.  
  543. /*function:external*/
  544. extern Ppoint3 ptk_subv3(C(Ppoint3 *) p1, C(Ppoint3 *) p2)
  545. PreANSI(Ppoint3 *p1)
  546. PreANSI(Ppoint3 *p2)
  547. /*
  548. ** \parambegin
  549. ** \param{Ppoint3 *}{p1}{3D vector}{IN}
  550. ** \param{Ppoint3 *}{p2}{3D vector}{IN}
  551. ** \paramend
  552. ** \blurb{This function evaluates the 3D vector \pardesc {v1-v2}.}
  553. */
  554. {
  555.   Ppoint3 v;
  556.  
  557.   v = ptk_point3(p1->x - p2->x, p1->y - p2->y, p1->z - p2->z);
  558.   return v;
  559. }  /* ptk_subv3 */
  560.  
  561. /*--------------------------------------------------------------------------*/
  562.  
  563. /*function:external*/
  564. extern Ppoint ptk_subv(C(Ppoint *) p1, C(Ppoint *) p2)
  565. PreANSI(Ppoint *p1)
  566. PreANSI(Ppoint *p2)
  567. /*
  568. ** \parambegin
  569. ** \param{Ppoint *}{p1}{2D vector}{IN}
  570. ** \param{Ppoint *}{p2}{2D vector}{IN}
  571. ** \paramend
  572. ** \blurb{This function evaluates the 2D vector \pardesc {v1-v2}.}
  573. */
  574. {
  575.   Ppoint v;
  576.  
  577.   v = ptk_point(p1->x - p2->x, p1->y - p2->y);
  578.   return v;
  579. }  /* ptk_subv */
  580.  
  581. /*--------------------------------------------------------------------------*/
  582.  
  583. /*function:external*/
  584. extern Ppoint3 ptk_addv3(C(Ppoint3 *) p1, C(Ppoint3 *) p2)
  585. PreANSI(Ppoint3 *p1)
  586. PreANSI(Ppoint3 *p2)
  587. /*
  588. ** \parambegin
  589. ** \param{Ppoint3 *}{p1}{3D vector}{IN}
  590. ** \param{Ppoint3 *}{p2}{3D vector}{IN}
  591. ** \paramend
  592. ** \blurb{This function adds the two 3D vectors \pardesc{v1} and \pardesc{v2},
  593. ** and returns the result.}
  594. */
  595. {
  596.   Ppoint3 v;
  597.  
  598.   v = ptk_point3(p1->x + p2->x, p1->y + p2->y, p1->z + p2->z);
  599.   return v;
  600. }  /* ptk_addv3 */
  601.  
  602. /*--------------------------------------------------------------------------*/
  603.  
  604. /*function:external*/
  605. extern Ppoint ptk_addv(C(Ppoint *) p1, C(Ppoint *) p2)
  606. PreANSI(Ppoint *p1)
  607. PreANSI(Ppoint *p2)
  608. /*
  609. ** \parambegin
  610. ** \param{Ppoint *}{p1}{2D vector}{IN}
  611. ** \param{Ppoint *}{p2}{2D vector}{IN}
  612. ** \paramend
  613. ** \blurb{This function adds the two 2D vectors \pardesc{v1}
  614. ** and \pardesc{v2}, and returns the result.}
  615. */
  616. {
  617.   Ppoint v;
  618.  
  619.   v = ptk_point(p1->x + p2->x, p1->y + p2->y);
  620.   return v;
  621. }  /* ptk_addv */
  622.  
  623. /*--------------------------------------------------------------------------*/
  624.  
  625. /*function:external*/
  626. extern void ptk_unitmatrix(C(Pmatrix) matrix)
  627. PreANSI(Pmatrix matrix)
  628. /*
  629. ** \parambegin
  630. ** \param{Pmatrix}{matrix}{3x3 matrix}{IN}
  631. ** \paramend
  632. ** \blurb{This procedure creates a unit $3\times 3$ matrix, and stores it in 
  633. ** \pardesc{matrix}.}
  634. */
  635. {
  636.   /*  */
  637.   Pint jj, ii;
  638.  
  639.   for (ii = 0; ii <= 2; ii++) 
  640.   {
  641.     for (jj = 0; jj <= 2; jj++)
  642.       matrix[ii][jj] = 0.0;
  643.     matrix[ii][ii] = 1.0;
  644.   }
  645. }  /* ptk_unitmatrix */
  646.  
  647. /*--------------------------------------------------------------------------*/
  648.  
  649. /*function:external*/
  650. extern void ptk_unitmatrix3(C(Pmatrix3) matrix)
  651. PreANSI(Pmatrix3 matrix)
  652. /*
  653. ** \parambegin
  654. ** \param{Pmatrix3}{matrix}{4x4 matrix}{IN}
  655. ** \paramend
  656. ** \blurb{This procedure creates a unit $4\times 4$ matrix, and stores it in 
  657. ** \pardesc{matrix}.}
  658. */
  659. {
  660.   /*  */
  661.   Pint jj, ii;
  662.  
  663.   for (ii = 0; ii <= 3; ii++) 
  664.   {
  665.     for (jj = 0; jj <= 3; jj++)
  666.       matrix[ii][jj] = 0.0;
  667.     matrix[ii][ii] = 1.0;
  668.   }
  669. }  /* ptk_unitmatrix3 */
  670.  
  671. /*--------------------------------------------------------------------------*/
  672.  
  673. /*function:external*/
  674. extern void ptk_transposematrix3(C(Pmatrix3) matrix, C(Pmatrix3) result)
  675. PreANSI(Pmatrix3 matrix)
  676. PreANSI(Pmatrix3 result)
  677. /*
  678. ** \parambegin
  679. ** \param{Pmatrix3}{matrix}{4x4 matrix}{IN}
  680. ** \param{Pmatrix3}{result}{4x4 matrix}{OUT}
  681. ** \paramend
  682. ** \blurb{This function transposes \pardesc{matrix}, and returns the resul
  683. ** in \pardesc{result}.
  684. ** Note that \pardesc{restut} can be the same variable
  685. ** as \pardesc{matrix} since a copy is made 
  686. ** first.}
  687. */
  688. {
  689.   Pmatrix3 temp;
  690.   Pint i, j;
  691.  
  692.   memcpy(temp, matrix, sizeof(Pmatrix3));
  693.   for (i = 0; i <= 3; i++) 
  694.   {
  695.     for (j = 0; j <= 3; j++)
  696.       result[i][j] = temp[j][i];
  697.   }
  698. }  /* ptk_transposematrix3 */
  699.  
  700. /*--------------------------------------------------------------------------*/
  701.  
  702. /*function:external*/
  703. extern void ptk_transposematrix(C(Pmatrix) matrix, C(Pmatrix) result)
  704. PreANSI(Pmatrix matrix)
  705. PreANSI(Pmatrix result)
  706. /*
  707. ** \parambegin
  708. ** \param{Pmatrix}{matrix}{3x3 matrix}{IN}
  709. ** \param{Pmatrix}{result}{3x3 matrix}{OUT}
  710. ** \paramend
  711. ** \blurb{This function transposes \pardesc{matrix}, and returns the resul
  712. ** in \pardesc{result}.
  713. ** Note that \pardesc{restut} can be the same variable
  714. ** as \pardesc{matrix} since a copy is made 
  715. ** first.}
  716. */
  717. {
  718.   Pmatrix temp;
  719.   Pint i, j;
  720.  
  721.   memcpy(temp, matrix, sizeof(Pmatrix));
  722.   for (i = 0; i <= 2; i++) 
  723.   {
  724.     for (j = 0; j <= 2; j++)
  725.       result[i][j] = temp[j][i];
  726.   }
  727. }  /* ptk_transposematrix */
  728.  
  729. /*--------------------------------------------------------------------------*/
  730.  
  731. /*function:external*/
  732. extern void ptk_multiplymatrix3(C(Pmatrix3) matrix1, C(Pmatrix3) matrix2, 
  733.                           C(Pmatrix3) result)
  734. PreANSI(Pmatrix3 matrix1)
  735. PreANSI(Pmatrix3 matrix2)
  736. PreANSI(Pmatrix3 result)
  737. /*
  738. ** \parambegin
  739. ** \param{Pmatrix3}{matrix1}{4x4 matrix}{IN}
  740. ** \param{Pmatrix3}{matrix2}{4x4 matrix}{IN}
  741. ** \param{Pmatrix3}{result}{4x4 matrix}{OUT}
  742. ** \paramend
  743. ** \blurb{This function makes \pardesc{result} the product of
  744. ** the $4 \times 4$ matrices  \pardesc{matrix1} and
  745. ** \pardesc{matrix2}, with {\tt result $\leftarrow$ matrix1 * matrix2}.
  746. ** Note that \pardesc{result} can also be \pardesc{matrix1} or \pardesc{matrix2}
  747. ** since a copy is made.
  748. ** This function uses Winograd's method, see "Handbook of Algorithms
  749. ** and Data Structures", G.H. Gonnet, Addison Wesley.}
  750. */
  751. {
  752.   Pint ii, jj;
  753.   Pmatrix3 sum;
  754.   Pfloat temp;
  755.   Pfloat d[4], e[4];
  756.  
  757.   for (ii = 0; ii <= 3; ii++) 
  758.   {
  759.     d[ii] = (matrix1[ii][1] * matrix1[ii][0]) + 
  760.             (matrix1[ii][3] * matrix1[ii][2]);
  761.     e[ii] = (matrix2[0][ii] * matrix2[1][ii]) + 
  762.             (matrix2[2][ii] * matrix2[3][ii]);
  763.   }
  764.  
  765.   for (ii = 0; ii <= 3; ii++) 
  766.   {
  767.     for (jj = 0; jj <= 3; jj++) 
  768.     {
  769.       temp = (matrix1[ii][1] + matrix2[0][jj]) * 
  770.              (matrix1[ii][0] + matrix2[1][jj]);
  771.       temp += (matrix1[ii][3] + matrix2[2][jj]) * 
  772.               (matrix1[ii][2] + matrix2[3][jj]) - d[ii] - e[jj];
  773.       sum[ii][jj] = temp;
  774.     }
  775.   }
  776.   /* copy matrix back in case result is one of the other two */
  777.   memcpy(result, sum, sizeof(Pmatrix3));
  778. }  /* ptk_multiplymatrix3 */
  779.  
  780. /*--------------------------------------------------------------------------*/
  781.  
  782. /*function:external*/
  783. extern void ptk_multiplymatrix(C(Pmatrix) matrix1, C(Pmatrix) matrix2, 
  784.                           C(Pmatrix) result)
  785. PreANSI(Pmatrix matrix1)
  786. PreANSI(Pmatrix matrix2)
  787. PreANSI(Pmatrix result)
  788. /*
  789. ** \parambegin
  790. ** \param{Pmatrix}{matrix1}{3x3 matrix}{IN}
  791. ** \param{Pmatrix}{matrix2}{3x3 matrix}{IN}
  792. ** \param{Pmatrix}{result}{3x3 matrix}{OUT}
  793. ** \paramend
  794. ** \blurb{This function makes \pardesc{result} the product of the
  795. ** $3 \times 3$ matrices \pardesc{matrix1} and
  796. ** \pardesc{matrix2}, with {\tt result $\leftarrow$ matrix1 * matrix2}.
  797. ** Note that \pardesc{result} can also be \pardesc{matrix1} or \pardesc{matrix2}
  798. ** since a copy is made.}
  799. */
  800. {
  801.   Pint ii, jj, kk;
  802.   Pmatrix sum;
  803.   Pfloat temp;
  804.   Pfloat d[3], e[3];
  805.  
  806.   for (ii = 0; ii <= 2; ii++)
  807.     for (jj = 0; jj <= 2; jj++)
  808.       sum[ii][jj] = 0.0;
  809.   for (ii = 0; ii <= 2; ii++)
  810.     for (jj = 0; jj <= 2; jj++)
  811.       for (kk = 0; kk <= 2; kk++)
  812.         sum[ii][jj] += (matrix1[ii][kk] * matrix2[kk][jj]);
  813.  
  814. /*  for (ii = 0; ii <= 2; ii++) 
  815.   {
  816.     d[ii] = (matrix1[ii][1] * matrix1[ii][0]);
  817.     e[ii] = (matrix2[0][ii] * matrix2[1][ii]); 
  818.   }
  819.   for (ii = 0; ii <= 2; ii++) 
  820.   {
  821.     for (jj = 0; jj <= 2; jj++) 
  822.     {
  823.       temp = (matrix1[ii][1] + matrix2[0][jj]) * 
  824.              (matrix1[ii][0] + matrix2[1][jj]);
  825.       temp += (matrix1[ii][2] * matrix2[2][2]) - d[ii] - e[jj];
  826.       sum[ii][jj] = temp;
  827.     }
  828.   }
  829. */
  830.   /* copy matrix back in case result is one of the other two */
  831.   memcpy(result, sum, sizeof(Pmatrix));
  832. }  /* ptk_multiplymatrix */
  833.  
  834. /*--------------------------------------------------------------------------*/
  835.  
  836. /*function:external*/
  837. extern void ptk_concatenatematrix3(C(Pcompose_type) operation, 
  838.             C(Pmatrix3) matrix1, C(Pmatrix3) matrix2, C(Pmatrix3) result)
  839. PreANSI(Pcompose_type operation)
  840. PreANSI(Pmatrix3 matrix1)
  841. PreANSI(Pmatrix3 matrix2)
  842. PreANSI(Pmatrix3 result)
  843. /*
  844. ** \parambegin
  845. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  846. ** \param{Pmatrix3}{matrix1}{4x4 matrix}{IN}
  847. ** \param{Pmatrix3}{matrix2}{4x4 matrix}{IN}
  848. ** \param{Pmatrix3}{result}{4x4 matrix}{OUT}
  849. ** \paramend
  850. ** \blurb{This function concatenates the $4 \times 4$ matrices
  851. ** \pardesc{matrix1}  and \pardesc{matrix2}
  852. ** on the basis of \pardesc{operation}.
  853. ** The result is stored in \pardesc{result}.
  854. ** Note that \pardesc{result} can also be \pardesc{matrix1} or \pardesc{matrix2}
  855. ** since a copy is made. When \pardesc{operation} is 
  856. ** \pardesc{preconcatenate}, then \pardesc{result $\leftarrow$ matrix2 * matrix1}.
  857. ** When \pardesc{operation} is \pardesc{postconcatenate}, 
  858. ** \pardesc{result $\leftarrow$  matrix1 * matrix2}.}
  859. */
  860. {
  861.   switch (operation) 
  862.   {
  863.   case PTYPE_PRECONCAT:  
  864.     ptk_multiplymatrix3(matrix2, matrix1, result);
  865.     break;
  866.  
  867.   case PTYPE_POSTCONCAT:  
  868.     ptk_multiplymatrix3(matrix1, matrix2, result);
  869.     break;
  870.  
  871.   case PTYPE_REPLACE:  
  872.     memcpy(result, matrix1, sizeof(Pmatrix3));
  873.     break;
  874.   }
  875. }  /* ptk_concatenatematrix3 */
  876.  
  877. /*--------------------------------------------------------------------------*/
  878.  
  879. /*function:external*/
  880. extern void ptk_concatenatematrix(C(Pcompose_type) operation, 
  881.             C(Pmatrix) matrix1, C(Pmatrix) matrix2, C(Pmatrix) result)
  882. PreANSI(Pcompose_type operation)
  883. PreANSI(Pmatrix matrix1)
  884. PreANSI(Pmatrix matrix2)
  885. PreANSI(Pmatrix result)
  886. /*
  887. ** \parambegin
  888. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  889. ** \param{Pmatrix}{matrix1}{3x3 matrix}{IN}
  890. ** \param{Pmatrix}{matrix2}{3x3 matrix}{IN}
  891. ** \param{Pmatrix}{result}{3x3 matrix}{OUT}
  892. ** \paramend
  893. ** \blurb{This function concatenates
  894. ** the $3 \times 3$ matrices \pardesc{matrix1}  and \pardesc{matrix2}
  895. ** on the basis of \pardesc{operation}.
  896. ** The result is stored in \pardesc{result}.
  897. ** Note that \pardesc{result} can also be \pardesc{matrix1} or \pardesc{matrix2}
  898. ** since a copy is made. When \pardesc{operation} is
  899. ** \pardesc{preconcatenate}, then \pardesc{result $\leftarrow$ matrix2 * matrix1}.
  900. ** When \pardesc{operation} is \pardesc{postconcatenate}, 
  901. ** \pardesc{result $\leftarrow$ matrix1 * matrix2}.}
  902. */
  903. {
  904.   switch (operation) 
  905.   {
  906.   case PTYPE_PRECONCAT:  
  907.     ptk_multiplymatrix(matrix2, matrix1, result);
  908.     break;
  909.  
  910.   case PTYPE_POSTCONCAT:  
  911.     ptk_multiplymatrix(matrix1, matrix2, result);
  912.     break;
  913.  
  914.   case PTYPE_REPLACE:  
  915.     memcpy(result, matrix1, sizeof(Pmatrix));
  916.     break;
  917.   }
  918. }  /* ptk_concatenatematrix */
  919.  
  920. /*--------------------------------------------------------------------------*/
  921.  
  922. /*function:external*/
  923. extern void ptk_shift3(C(Ppoint3 *) shift, C(Pcompose_type) operation, 
  924.                        C(Pmatrix3) matrix)
  925. PreANSI(Ppoint3 *shift)
  926. PreANSI(Pcompose_type operation)
  927. PreANSI(Pmatrix3 matrix)
  928. /*
  929. ** \parambegin
  930. ** \param{Ppoint3 *}{shift}{shift factor}{IN}
  931. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  932. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  933. ** \paramend
  934. ** \blurb{This function computes
  935. **  a matrix to perform the specified 3D shift and concatenates 
  936. ** this matrix with \pardesc{matrix} on the basis of \pardesc{operation}.}
  937. */
  938. {
  939.   Pmatrix3 temp;
  940.  
  941.   ptk_unitmatrix3(temp);
  942.   temp[0][3] = shift->x;
  943.   temp[1][3] = shift->y;
  944.   temp[2][3] = shift->z;
  945.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  946. }  /* ptk_shift3 */
  947.  
  948. /*--------------------------------------------------------------------------*/
  949.  
  950. /*function:external*/
  951. extern void ptk_shift(C(Ppoint *) shift, C(Pcompose_type) operation, 
  952.                        C(Pmatrix) matrix)
  953. PreANSI(Ppoint *shift)
  954. PreANSI(Pcompose_type operation)
  955. PreANSI(Pmatrix matrix)
  956. /*
  957. ** \parambegin
  958. ** \param{Ppoint *}{shift}{shift factor}{IN}
  959. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  960. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  961. ** \paramend
  962. ** \blurb{This function computes
  963. ** a matrix to perform the specified 2D shift and concatenates 
  964. ** this matrix with \pardesc{matrix} on the basis of \pardesc{operation}.}
  965. */
  966. {
  967.   Pmatrix temp;
  968.  
  969.   ptk_unitmatrix(temp);
  970.   temp[0][2] = shift->x;
  971.   temp[1][2] = shift->y;
  972.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  973. }  /* ptk_shift */
  974.  
  975. /*--------------------------------------------------------------------------*/
  976.  
  977. /*function:external*/
  978. extern void ptk_scale3(C(Ppoint3 *) scale, C(Pcompose_type) operation, 
  979.                        C(Pmatrix3) matrix)
  980. PreANSI(Ppoint3 *scale)
  981. PreANSI(Pcompose_type operation)
  982. PreANSI(Pmatrix3 matrix)
  983. /*
  984. ** \parambegin
  985. ** \param{Ppoint3 *}{shift}{shift factor}{IN}
  986. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  987. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  988. ** \paramend
  989. ** \blurb{This function 
  990. ** computes a matrix to perform the specified 3D scale and concatenates 
  991. ** this with \pardesc{matrix} on the basis of \pardesc{operation}.}
  992. */
  993. {
  994.   Pmatrix3 temp;
  995.  
  996.   ptk_unitmatrix3(temp);
  997.   temp[0][0] = scale->x;  
  998.   temp[1][1] = scale->y;
  999.   temp[2][2] = scale->z;
  1000.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  1001. }  /* ptk_scale3 */
  1002.  
  1003. /*--------------------------------------------------------------------------*/
  1004.  
  1005. /*function:external*/
  1006. extern void ptk_scale(C(Ppoint *) scale, C(Pcompose_type) operation, 
  1007.                        C(Pmatrix) matrix)
  1008. PreANSI(Ppoint *scale)
  1009. PreANSI(Pcompose_type operation)
  1010. PreANSI(Pmatrix matrix)
  1011. /*
  1012. ** \parambegin
  1013. ** \param{Ppoint *}{shift}{shift factor}{IN}
  1014. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1015. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  1016. ** \paramend
  1017. ** \blurb{This function computes
  1018. **  a matrix to perform the specified 2D scale and concatenates 
  1019. ** this matrix with \pardesc{matrix} on the basis of \pardesc{operation}.}
  1020. */
  1021. {
  1022.   Pmatrix temp;
  1023.  
  1024.   ptk_unitmatrix(temp);
  1025.   temp[0][0] = scale->x;  
  1026.   temp[1][1] = scale->y;
  1027.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  1028. }  /* ptk_scale */
  1029.  
  1030. /*--------------------------------------------------------------------------*/
  1031.  
  1032. /*function:external*/
  1033. extern void ptk_rotatecs3(C(Pfloat) costheta, C(Pfloat) sinetheta, 
  1034.                           C(ptkeaxistype) axis, C(Pcompose_type) operation, 
  1035.                           C(Pmatrix3) matrix)
  1036. PreANSI(Pfloat costheta)
  1037. PreANSI(Pfloat sinetheta)
  1038. PreANSI(ptkeaxistype axis)
  1039. PreANSI(Pcompose_type operation)
  1040. PreANSI(Pmatrix3 matrix)
  1041. /*
  1042. ** \parambegin
  1043. ** \param{Pfloat}{costheta}{cosine of angle}{IN}
  1044. ** \param{Pfloat}{sinetheta}{sine of angle}{IN}
  1045. ** \param{ptkeaxistype}{axis}{x, y or z axis}{IN}
  1046. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1047. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  1048. ** \paramend
  1049. ** \blurb{This function computes
  1050. ** a matrix to perform the specified 3D rotation and concatenates 
  1051. ** this matrix with 
  1052. ** \pardesc{matrix} on the basis of \pardesc{operation}.
  1053. ** This form assumes that the rotation is specified
  1054. ** using the $\cos(theta)$ and $\sin(theta)$ terms.
  1055. ** Note that no check is made to ensure that the sum of the squares of these
  1056. ** terms is 1.}
  1057. */
  1058. {
  1059.   Pmatrix3 temp;
  1060.  
  1061.   ptk_unitmatrix3(temp);
  1062.   switch (axis) 
  1063.   {
  1064.   case PTKEXAXIS:
  1065.     temp[2][2] = costheta;
  1066.     temp[1][1] = costheta;
  1067.     temp[1][2] = -sinetheta;
  1068.     temp[2][1] = sinetheta;
  1069.     break;
  1070.  
  1071.   case PTKEYAXIS:
  1072.     temp[0][0] = costheta;
  1073.     temp[2][2] = costheta;
  1074.     temp[0][2] = sinetheta;
  1075.     temp[2][0] = -sinetheta;
  1076.     break;
  1077.  
  1078.   case PTKEZAXIS:
  1079.     temp[0][0] = costheta;
  1080.     temp[1][1] = costheta;
  1081.     temp[0][1] = -sinetheta;
  1082.     temp[1][0] = sinetheta;
  1083.     break;
  1084.   }
  1085.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  1086. }  /* ptk_rotatecs3 */
  1087.  
  1088. /*--------------------------------------------------------------------------*/
  1089.  
  1090. /*function:external*/
  1091. extern void ptk_rotatecs(C(Pfloat) costheta, C(Pfloat) sinetheta, 
  1092.                           C(ptkeaxistype) axis, C(Pcompose_type) operation, 
  1093.                           C(Pmatrix) matrix)
  1094. PreANSI(Pfloat costheta)
  1095. PreANSI(Pfloat sinetheta)
  1096. PreANSI(ptkeaxistype axis)
  1097. PreANSI(Pcompose_type operation)
  1098. PreANSI(Pmatrix matrix)
  1099. /*
  1100. ** \parambegin
  1101. ** \param{Pfloat}{costheta}{cosine of angle}{IN}
  1102. ** \param{Pfloat}{sinetheta}{sine of angle}{IN}
  1103. ** \param{ptkeaxistype}{axis}{x, y or z axis}{IN}
  1104. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1105. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  1106. ** \paramend
  1107. ** \blurb{This function computes
  1108. ** a matrix to perform the specified 2D rotation and concatenates 
  1109. ** this matrix with 
  1110. ** \pardesc{matrix} on the basis of \pardesc{operation}.
  1111. ** This form assumes that the rotation is specified
  1112. ** using the $\cos(theta)$ and $\sin(theta)$ terms.
  1113. ** Note that no check is made to ensure that the sum of the squares of these
  1114. ** terms is 1.}
  1115. */
  1116. {
  1117.   Pmatrix temp;
  1118.  
  1119.   ptk_unitmatrix(temp);
  1120.   switch (axis) 
  1121.   {
  1122.   case PTKEXAXIS:
  1123.     temp[2][2] = costheta;
  1124.     temp[1][1] = costheta;
  1125.     temp[1][2] = -sinetheta;
  1126.     temp[2][1] = sinetheta;
  1127.     break;
  1128.  
  1129.   case PTKEYAXIS:
  1130.     temp[0][0] = costheta;
  1131.     temp[2][2] = costheta;
  1132.     temp[0][2] = sinetheta;
  1133.     temp[2][0] = -sinetheta;
  1134.     break;
  1135.  
  1136.   case PTKEZAXIS:
  1137.     temp[0][0] = costheta;
  1138.     temp[1][1] = costheta;
  1139.     temp[0][1] = -sinetheta;
  1140.     temp[1][0] = sinetheta;
  1141.     break;
  1142.   }
  1143.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  1144. }  /* ptk_rotatecs */
  1145.  
  1146. /*--------------------------------------------------------------------------*/
  1147.  
  1148. /*function:external*/
  1149. extern void ptk_rotate3(C(Pfloat) rotation, C(ptkeaxistype) axis, 
  1150.                         C(Pcompose_type) operation, C(Pmatrix3) matrix)
  1151. PreANSI(Pfloat rotation)
  1152. PreANSI(ptkeaxistype axis)
  1153. PreANSI(Pcompose_type operation) 
  1154. PreANSI(Pmatrix3 matrix)
  1155. /*
  1156. ** \parambegin
  1157. ** \param{Pfloat}{rotation}{angle in degrees}{IN}
  1158. ** \param{ptkeaxistype}{axis}{x, y or z axis}{IN}
  1159. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1160. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  1161. ** \paramend
  1162. ** \blurb{This function computes
  1163. ** a matrix to perform the specified 3D rotation and concatenates 
  1164. ** this matrix with 
  1165. ** \pardesc{matrix} on the basis of \pardesc{operation}.
  1166. ** \pardesc{rotation} is expressed in degrees.}
  1167. */
  1168. {
  1169.   Pfloat costheta;
  1170.   Pfloat sinetheta;
  1171.  
  1172.   costheta  = cos((rotation * ptkcdtor));
  1173.   sinetheta = sin((rotation * ptkcdtor));
  1174.   ptk_rotatecs3(costheta, sinetheta, axis, operation, matrix);
  1175. }  /* ptk_rotate3 */
  1176.  
  1177. /*--------------------------------------------------------------------------*/
  1178.  
  1179. /*function:external*/
  1180. extern void ptk_rotate(C(Pfloat) rotation, C(ptkeaxistype) axis, 
  1181.                         C(Pcompose_type) operation, C(Pmatrix) matrix)
  1182. PreANSI(Pfloat rotation)
  1183. PreANSI(ptkeaxistype axis)
  1184. PreANSI(Pcompose_type operation) 
  1185. PreANSI(Pmatrix matrix)
  1186. /*
  1187. ** \parambegin
  1188. ** \param{Pfloat}{rotation}{angle in degrees}{IN}
  1189. ** \param{ptkeaxistype}{axis}{x, y or z axis}{IN}
  1190. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1191. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  1192. ** \paramend
  1193. ** \blurb{This function computes
  1194. ** a matrix to perform the specified 2D rotation and concatenates 
  1195. ** this matrix with 
  1196. ** \pardesc{matrix} on the basis of \pardesc{operation}.
  1197. ** \pardesc{rotation} is expressed in degrees.}
  1198. */
  1199. {
  1200.   Pfloat costheta;
  1201.   Pfloat sinetheta;
  1202.  
  1203.   costheta  = cos((rotation * ptkcdtor));
  1204.   sinetheta = sin((rotation * ptkcdtor));
  1205.   ptk_rotatecs(costheta, sinetheta, axis, operation, matrix);
  1206. }  /* ptk_rotate */
  1207.  
  1208. /*--------------------------------------------------------------------------*/
  1209.  
  1210. /*function:external*/
  1211. extern void ptk_shear3(C(ptkeaxistype) shearaxis, C(ptkeaxistype) sheardir, 
  1212.                        C(Pfloat) shearfactor, C(Pcompose_type) operation, 
  1213.                        C(Pmatrix3) matrix)
  1214. PreANSI(ptkeaxistype shearaxis)
  1215. PreANSI(ptkeaxistype sheardir)
  1216. PreANSI(Pfloat shearfactor)
  1217. PreANSI(Pcompose_type operation)
  1218. PreANSI(Pmatrix3 matrix)
  1219. /*
  1220. ** \parambegin
  1221. ** \param{ptkeaxistype}{shearaxis}{x, y or z axis}{IN}
  1222. ** \param{ptkeaxistype}{sheardir}{x, y or z direction}{IN}
  1223. ** \param{Pfloat}{shearfactor}{shear factor}{IN}
  1224. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1225. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  1226. ** \paramend
  1227. ** \blurb{This function computes a matrix to perform the specified 3D
  1228. **  shear and concatenates 
  1229. ** this matrix with 
  1230. ** \pardesc{matrix} on the basis of \pardesc{operation}.
  1231. ** The shear is specified as an amount \pardesc{f} about axis \pardesc{i}
  1232. ** in direction \pardesc{j}.}
  1233. */
  1234. {
  1235.   Pmatrix3 temp;
  1236.  
  1237.   ptk_unitmatrix3(temp);
  1238.   temp[(Pint)shearaxis - 1][(Pint)sheardir - 1] = shearfactor;
  1239.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  1240. }  /* ptk_shear3 */
  1241.  
  1242. /*--------------------------------------------------------------------------*/
  1243.  
  1244. /*function:external*/
  1245. extern void ptk_shear(C(ptkeaxistype) shearaxis, C(ptkeaxistype) sheardir, 
  1246.                        C(Pfloat) shearfactor, C(Pcompose_type) operation, 
  1247.                        C(Pmatrix) matrix)
  1248. PreANSI(ptkeaxistype shearaxis)
  1249. PreANSI(ptkeaxistype sheardir)
  1250. PreANSI(Pfloat shearfactor)
  1251. PreANSI(Pcompose_type operation)
  1252. PreANSI(Pmatrix matrix)
  1253. /*
  1254. ** \parambegin
  1255. ** \param{ptkeaxistype}{shearaxis}{x or y axis}{IN}
  1256. ** \param{ptkeaxistype}{sheardir}{x or y direction}{IN}
  1257. ** \param{Pfloat}{shearfactor}{shear factor}{IN}
  1258. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1259. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  1260. ** \paramend
  1261. ** \blurb{This function computes a matrix to perform the specified 2D
  1262. **  shear and concatenates 
  1263. ** this matrix with 
  1264. ** \pardesc{matrix} on the basis of \pardesc{operation}.
  1265. ** The shear is specified as an amount \pardesc{f} about axis \pardesc{i}
  1266. ** in direction \pardesc{j}.}
  1267. */
  1268. {
  1269.   Pmatrix temp;
  1270.  
  1271.   ptk_unitmatrix(temp);
  1272.   temp[(Pint)shearaxis - 1][(Pint)sheardir - 1] = shearfactor;
  1273.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  1274. }  /* ptk_shear */
  1275.  
  1276. /*--------------------------------------------------------------------------*/
  1277.  
  1278. /*function:external*/
  1279. extern void ptk_rotatevv3(C(Ppoint3 *) v1, C(Ppoint3 *) v2, 
  1280.                           C(Pcompose_type) operation, 
  1281.                           C(Pmatrix3) matrix, C(Pint *)error)
  1282. PreANSI(Ppoint3 *v1)
  1283. PreANSI(Ppoint3 *v2)
  1284. PreANSI(Pcompose_type operation) 
  1285. PreANSI(Pmatrix3 matrix)
  1286. PreANSI(Pint *error)
  1287. /*
  1288. ** \parambegin
  1289. ** \param{Ppoint3 *}{v1}{3D vector}{IN}
  1290. ** \param{Ppoint3 *}{v2}{3D vector}{IN}
  1291. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1292. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  1293. ** \param{Pint *}{error}{error code}{OUT}
  1294. ** \paramend
  1295. ** \blurb{This function computes
  1296. **  a matrix to perform the rotation (about the origin) of the 3D vector
  1297. ** \pardesc{v1} to the 3D vector 
  1298. ** \pardesc{v2}, and concatenates this matrix with 
  1299. ** \pardesc{matrix} on the basis of \pardesc{operation} (\cite{rog:mecg}, 
  1300. **  pages 55--59). If the parameters are invalid, \pardesc{error} is set to 
  1301. ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
  1302. */
  1303. {
  1304.   Pmatrix3 temp;
  1305.   Ppoint3 dv1, dv2, nvec;
  1306.   Pfloat t1, costheta, sinetheta;
  1307.  
  1308.   /* Create a rotation matrix (about origin) to rotate from
  1309.   ** vector dv1 to vector dv2
  1310.   ** Cf ROGERS and ADAMS page 55-59
  1311.   */
  1312.  
  1313.   *error = ptkcpcok;   /* Assume OK */
  1314.   dv1 = ptk_unitv3(v1);
  1315.   dv2 = ptk_unitv3(v2);
  1316.  
  1317.   /* Get sine theta (cross product) and cos theta (dot product) */
  1318.  
  1319.   costheta = dv1.x * dv2.x + dv1.y * dv2.y + dv1.z * dv2.z;
  1320.   nvec = ptk_crossv3(&dv1, &dv2);
  1321.  
  1322.   /* How do I choose sign of sinetheta?
  1323.   ** let A, B and K be unit vectors, then
  1324.   ** A X B :=  sin(theta) K
  1325.   ** The sign of sinetheta is not important
  1326.   ** because of the way it appears in the subsequent formulae
  1327.   */
  1328.  
  1329.   sinetheta = ptk_modv3(&nvec);
  1330.   if (ptk_equal(fabs(sinetheta), 0.0)) 
  1331.   {
  1332.     *error = -1;
  1333.     sinetheta = ptkcpceps;
  1334.   }
  1335.   nvec = ptk_unitv3(&nvec);
  1336.   ptk_unitmatrix3(temp);   /* get a unit matrix */
  1337.   t1 = 1.0 - costheta;
  1338.   temp[0][0] = (nvec.x * nvec.x) + (1.0 - (nvec.x * nvec.x)) * costheta;
  1339.   temp[1][1] = (nvec.y * nvec.y) + (1.0 - (nvec.y * nvec.y)) * costheta;
  1340.   temp[2][2] = (nvec.z * nvec.z) + (1.0 - (nvec.z * nvec.z)) * costheta;
  1341.   temp[1][0] = (nvec.x * nvec.y * t1) + (nvec.z * sinetheta);
  1342.   temp[2][0] = (nvec.x * nvec.z * t1) - (nvec.y * sinetheta); 
  1343.   temp[0][1] = (nvec.x * nvec.y * t1) - (nvec.z * sinetheta);
  1344.   temp[2][1] = (nvec.y * nvec.z * t1) + (nvec.x * sinetheta);
  1345.   temp[0][2] = (nvec.x * nvec.z * t1) + (nvec.y * sinetheta);
  1346.   temp[1][2] = (nvec.y * nvec.z * t1) - (nvec.x * sinetheta);
  1347.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  1348. }  /* ptk_rotatevv3 */
  1349.  
  1350. /*--------------------------------------------------------------------------*/
  1351.  
  1352. /*function:external*/
  1353. extern void ptk_rotatevv(C(Ppoint *) v1, C(Ppoint *) v2, 
  1354.                           C(Pcompose_type) operation, 
  1355.                           C(Pmatrix) matrix, C(Pint *)error)
  1356. PreANSI(Ppoint *v1)
  1357. PreANSI(Ppoint *v2)
  1358. PreANSI(Pcompose_type operation) 
  1359. PreANSI(Pmatrix matrix)
  1360. PreANSI(Pint *error)
  1361. /*
  1362. ** \parambegin
  1363. ** \param{Ppoint *}{v1}{2D vector}{IN}
  1364. ** \param{Ppoint *}{v2}{2D vector}{IN}
  1365. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1366. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  1367. ** \param{Pint *}{error}{error code}{OUT}
  1368. ** \paramend
  1369. ** \blurb{This function computes
  1370. **  a matrix to perform the rotation (about the origin) of the 3D vector
  1371. ** \pardesc{v1} to the 3D vector 
  1372. ** \pardesc{v2}, and concatenates this matrix with 
  1373. ** \pardesc{matrix} on the basis of \pardesc{operation} (\cite{rog:mecg}, 
  1374. **  pages 55--59). If the parameters are invalid, \pardesc{error} is set to 
  1375. ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
  1376. */
  1377. {
  1378.   Pmatrix temp;
  1379.   Ppoint dv1, dv2, nvec;
  1380.   Pfloat t1, costheta, sinetheta;
  1381.  
  1382.   /* Create a rotation matrix (about origin) to rotate from
  1383.   ** vector dv1 to vector dv2
  1384.   ** Cf ROGERS and ADAMS page 55-59
  1385.   */
  1386.  
  1387.   *error = ptkcpcok;   /* Assume OK */
  1388.   dv1 = ptk_unitv(v1);
  1389.   dv2 = ptk_unitv(v2);
  1390.  
  1391.   /* Get sine theta (cross product) and cos theta (dot product) */
  1392.  
  1393.   costheta = dv1.x * dv2.x + dv1.y * dv2.y;
  1394.   nvec = ptk_point(0.0, 0.0);
  1395.  
  1396.   /* How do I choose sign of sinetheta?
  1397.   ** let A, B and K be unit vectors, then
  1398.   ** A X B :=  sin(theta) K
  1399.   ** The sign of sinetheta is not important
  1400.   ** because of the way it appears in the subsequent formulae
  1401.   */
  1402.  
  1403.   sinetheta = ptk_modv(&nvec);
  1404.   if (ptk_equal(fabs(sinetheta), 0.0)) 
  1405.   {
  1406.     *error = -1;
  1407.     sinetheta = ptkcpceps;
  1408.   }
  1409.   nvec = ptk_unitv(&nvec);
  1410.   ptk_unitmatrix(temp);   /* get a unit matrix */
  1411.   t1 = 1.0 - costheta;
  1412.   temp[0][0] = (nvec.x * nvec.x) + (1.0 - (nvec.x * nvec.x)) * costheta;
  1413.   temp[1][1] = (nvec.y * nvec.y) + (1.0 - (nvec.y * nvec.y)) * costheta;
  1414.   temp[2][2] = costheta;
  1415.   temp[1][0] = (nvec.x * nvec.y * t1);
  1416.   temp[2][0] = (nvec.y * sinetheta); 
  1417.   temp[0][1] = (nvec.x * nvec.y * t1);
  1418.   temp[2][1] = (nvec.x * sinetheta);
  1419.   temp[0][2] = (nvec.y * sinetheta);
  1420.   temp[1][2] = (nvec.x * sinetheta);
  1421.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  1422. }  /* ptk_rotatevv */
  1423.  
  1424. /*--------------------------------------------------------------------------*/
  1425.  
  1426. /*function:external*/
  1427. extern void ptk_rotateline3(C(Ppoint3 *) p1, C(Ppoint3 *) p2, 
  1428.                             C(Pfloat) theta, C(Pcompose_type) operation, 
  1429.                             C(Pmatrix3) matrix, C(Pint *)error)
  1430. PreANSI(Ppoint3 *p1)
  1431. PreANSI(Ppoint3 *p2)
  1432. PreANSI(Pfloat theta) 
  1433. PreANSI(Pcompose_type operation)
  1434. PreANSI(Pmatrix3 matrix)
  1435. PreANSI(Pint *error)
  1436. /*
  1437. ** \parambegin
  1438. ** \param{Ppoint3 *}{p1}{3D point}{IN}
  1439. ** \param{Ppoint3 *}{p2}{3D point}{IN}
  1440. ** \param{Pfloat}{theta}{angle in degrees}{IN}
  1441. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1442. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  1443. ** \param{Pint *}{error}{error code}{OUT}
  1444. ** \paramend
  1445. ** \blurb{This function computes
  1446. ** a matrix to perform a 3D rotation of
  1447. ** \pardesc{theta} degrees
  1448. ** about the line connecting \pardesc{p1} to \pardesc{p2},
  1449. ** and concatenates this matrix with 
  1450. ** \pardesc{matrix} on the basis of \pardesc{operation} (\cite{rog:mecg}, 
  1451. **  pages 55--59).
  1452. ** If the parameters are invalid, \pardesc{error} is set to 
  1453. ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
  1454. */
  1455. {
  1456.   Pmatrix3 temp;
  1457.   Ppoint3 nvec, mnvec;
  1458.   Pfloat t1, costheta, sinetheta;
  1459.  
  1460.   *error = ptkcpcok;
  1461.   /* make a unit vector, i.e. direction cosines
  1462.   ** nvec = p2 - p1.
  1463.   */
  1464.   nvec = ptk_subv3(p2, p1);
  1465.   if (ptk_nullv3(&nvec)) {*error = -1;}
  1466.   nvec = ptk_unitv3(&nvec);
  1467.   costheta = cos(theta * ptkcdtor);
  1468.   sinetheta = sin(theta * ptkcdtor);
  1469.   t1 = 1.0 - costheta;
  1470.   ptk_unitmatrix3(temp); /* get a unit matrix */
  1471.   temp[0][0] = (nvec.x * nvec.x) + (1.0 - (nvec.x * nvec.x)) * costheta;
  1472.   temp[1][1] = (nvec.y * nvec.y) + (1.0 - (nvec.y * nvec.y)) * costheta;
  1473.   temp[2][2] = (nvec.z * nvec.z) + (1.0 - (nvec.z * nvec.z)) * costheta;
  1474.   temp[1][0] = (nvec.x * nvec.y * t1) + (nvec.z * sinetheta);
  1475.   temp[2][0] = (nvec.x * nvec.z * t1) - (nvec.y * sinetheta); 
  1476.   temp[0][1] = (nvec.x * nvec.y * t1) - (nvec.z * sinetheta);
  1477.   temp[2][1] = (nvec.y * nvec.z * t1) + (nvec.x * sinetheta);
  1478.   temp[0][2] = (nvec.x * nvec.z * t1) + (nvec.y * sinetheta);
  1479.   temp[1][2] = (nvec.y * nvec.z * t1) - (nvec.x * sinetheta);
  1480.   mnvec = ptk_scalev3(p1, -1.0);
  1481.   ptk_shift3(&mnvec, PTYPE_PRECONCAT, temp);
  1482.   ptk_shift3(p1, PTYPE_POSTCONCAT, temp);
  1483.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  1484. }  /* ptk_rotateline3 */
  1485.  
  1486. /*--------------------------------------------------------------------------*/
  1487.  
  1488. /*function:external*/
  1489. extern void ptk_rotateline(C(Ppoint *) p1, C(Ppoint *) p2, 
  1490.                             C(Pfloat) theta, C(Pcompose_type) operation, 
  1491.                             C(Pmatrix) matrix, C(Pint *)error)
  1492. PreANSI(Ppoint *p1)
  1493. PreANSI(Ppoint *p2)
  1494. PreANSI(Pfloat theta) 
  1495. PreANSI(Pcompose_type operation)
  1496. PreANSI(Pmatrix matrix)
  1497. PreANSI(Pint *error)
  1498. /*
  1499. ** \parambegin
  1500. ** \param{Ppoint *}{p1}{2D point}{IN}
  1501. ** \param{Ppoint *}{p2}{2D point}{IN}
  1502. ** \param{Pfloat}{theta}{angle in degrees}{IN}
  1503. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1504. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  1505. ** \param{Pint *}{error}{error code}{OUT}
  1506. ** \paramend
  1507. ** \blurb{This function computes
  1508. ** a matrix to perform a 2D rotation of
  1509. ** \pardesc{theta} degrees
  1510. ** about the line connecting \pardesc{p1} to \pardesc{p2},
  1511. ** and concatenates this matrix with 
  1512. ** \pardesc{matrix} on the basis of \pardesc{operation}.
  1513. ** If the parameters are invalid, \pardesc{error} is set to 
  1514. ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
  1515. */
  1516. {
  1517.   Pmatrix temp;
  1518.   Ppoint nvec, mnvec;
  1519.   Pfloat t1, costheta, sinetheta;
  1520.  
  1521.   *error = ptkcpcok;
  1522.   /* make a unit vector, i.e. direction cosines
  1523.   ** nvec = p2 - p1.
  1524.   */
  1525.   nvec = ptk_subv(p2, p1);
  1526.   if (ptk_nullv(&nvec)) {*error = -1;}
  1527.   nvec = ptk_unitv(&nvec);
  1528.   costheta = cos(theta * ptkcdtor);
  1529.   sinetheta = sin(theta * ptkcdtor);
  1530.   t1 = 1.0 - costheta;
  1531.   ptk_unitmatrix(temp); /* get a unit matrix */
  1532.   temp[0][0] = (nvec.x * nvec.x) + (1.0 - (nvec.x * nvec.x)) * costheta;
  1533.   temp[1][1] = (nvec.y * nvec.y) + (1.0 - (nvec.y * nvec.y)) * costheta;
  1534.   temp[2][2] = costheta;
  1535.   temp[1][0] = (nvec.x * nvec.y * t1);
  1536.   temp[2][0] = (nvec.y * sinetheta); 
  1537.   temp[0][1] = (nvec.x * nvec.y * t1);
  1538.   temp[2][1] = (nvec.x * sinetheta);
  1539.   temp[0][2] = (nvec.y * sinetheta);
  1540.   temp[1][2] = (nvec.x * sinetheta);
  1541.   mnvec = ptk_scalev(p1, -1.0);
  1542.   ptk_shift(&mnvec, PTYPE_PRECONCAT, temp);
  1543.   ptk_shift(p1, PTYPE_POSTCONCAT, temp);
  1544.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  1545. }  /* ptk_rotateline */
  1546.  
  1547. /*--------------------------------------------------------------------------*/
  1548.  
  1549. /*function:external*/
  1550. extern Ppoint4 ptk_pt3topt4(C(Ppoint3 *) point)
  1551. PreANSI(Ppoint3 *point)
  1552. /*
  1553. ** \parambegin
  1554. ** \param{Ppoint3 *}{point}{3D point}{IN}
  1555. ** \paramend
  1556. ** \blurb{This function converts the 3D point \pardesc{point} to a 4D point,
  1557. ** which is returned as the result of the function. The $w$ coordinate of the 
  1558. ** 4D point is set to $1.0$.}
  1559. */
  1560. {
  1561.   Ppoint4 temp;
  1562.  
  1563.   temp.x = point->x;
  1564.   temp.y = point->y;
  1565.   temp.z = point->z;
  1566.   temp.w = 1.0;
  1567.   return temp;
  1568. }  /* ptk_pt3topt4 */
  1569.  
  1570. /*--------------------------------------------------------------------------*/
  1571.  
  1572. /*function:external*/
  1573. extern Ppoint3 ptk_pt4topt3(C(Ppoint4 *) point)
  1574. PreANSI(Ppoint4 *point)
  1575. /*
  1576. ** \parambegin
  1577. ** \param{Ppoint4 *}{point}{4D point}{IN}
  1578. ** \paramend
  1579. ** \blurb{This function converts the 4D point \pardesc{point} to a 3D point, by
  1580. ** dividing by $w$. The 3D point is returned as the result of the function.}
  1581. */
  1582. {
  1583.   Ppoint3 temp;
  1584.   Pfloat w;
  1585.  
  1586.   if (ptk_equal(fabs(point->w), 0.0))
  1587.     w = 1.0 / ptkcpceps;
  1588.   else
  1589.     w = 1.0 / point->w;
  1590.   temp.x = point->x * w;
  1591.   temp.y = point->y * w;
  1592.   temp.z = point->z * w;
  1593.   return temp;
  1594. }  /* ptk_pt4topt3 */
  1595.  
  1596. /*--------------------------------------------------------------------------*/
  1597.  
  1598. /*function:external*/
  1599. extern Ppoint4 ptk_transform4(C(Pmatrix3) matrix, C(Ppoint4 *) point)
  1600. PreANSI(Pmatrix3 matrix)
  1601. PreANSI(Ppoint4 *point)
  1602. /*
  1603. ** \parambegin
  1604. ** \param{Pmatrix3}{matrix}{4x4 matrix}{IN}
  1605. ** \param{Ppoint4 *}{point}{4D point}{IN}
  1606. ** \paramend
  1607. ** \blurb{This function performs the 4D transformation 
  1608. ** (that is, with no homogeneous division) of the
  1609. ** point \pardesc{pt} by the $4 \times 4$ matrix
  1610. ** \pardesc{matrix}, and returns the result as
  1611. ** the value of the function.}
  1612. */
  1613. {
  1614.   Ppoint4 temp;
  1615.  
  1616.   temp.x = matrix[0][0] * point->x + matrix[0][1] * point->y +
  1617.          matrix[0][2] * point->z + matrix[0][3] * point->w;
  1618.   temp.y = matrix[1][0] * point->x + matrix[1][1] * point->y +
  1619.          matrix[1][2] * point->z + matrix[1][3] * point->w;
  1620.   temp.z = matrix[2][0] * point->x + matrix[2][1] * point->y +
  1621.          matrix[2][2] * point->z + matrix[2][3] * point->w;
  1622.   temp.w = matrix[3][0] * point->x + matrix[3][1] * point->y +
  1623.          matrix[3][2] * point->z + matrix[3][3] * point->w;
  1624.   return temp;
  1625. }  /* ptk_transform4 */
  1626.  
  1627. /*--------------------------------------------------------------------------*/
  1628.  
  1629. /*function:external*/
  1630. extern Ppoint3 ptk_transform3(C(Pmatrix3) matrix, C(Ppoint3 *) point)
  1631. PreANSI(Pmatrix3 matrix)
  1632. PreANSI(Ppoint3 *point)
  1633. /*
  1634. ** \parambegin
  1635. ** \param{Pmatrix3}{matrix}{4x4 matrix}{IN}
  1636. ** \param{Ppoint3 *}{point}{3D point}{IN}
  1637. ** \paramend
  1638. ** \blurb{This function transforms the 3D point \pardesc{point} by
  1639. ** the $4 \times 4$ matrix \pardesc{matrix},
  1640. ** including homogeneous division.}
  1641. */
  1642. {
  1643.   Ppoint3 temp;
  1644.   Pfloat w;
  1645.  
  1646.   temp.x = matrix[0][0] * point->x + matrix[0][1] * point->y +
  1647.          matrix[0][2] * point->z + matrix[0][3];
  1648.   temp.y = matrix[1][0] * point->x + matrix[1][1] * point->y +
  1649.          matrix[1][2] * point->z + matrix[1][3];
  1650.   temp.z = matrix[2][0] * point->x + matrix[2][1] * point->y +
  1651.          matrix[2][2] * point->z + matrix[2][3];
  1652.   w = matrix[3][0] * point->x + matrix[3][1] * point->y + 
  1653.         matrix[3][2] * point->z + matrix[3][3];
  1654.   if (ptk_equal(fabs(w), 0.0))
  1655.     w = 1.0 / ptkcpceps;
  1656.   else
  1657.     w = 1.0 / w;
  1658.   temp.x *= w;
  1659.   temp.y *= w;
  1660.   temp.z *= w;
  1661.   return temp;
  1662. }  /* ptk_transform3 */
  1663.  
  1664. /*-------------------------------------------------------------------*/
  1665.  
  1666. /*function:external*/
  1667. extern Ppoint ptk_transform(C(Pmatrix) matrix, C(Ppoint *) point)
  1668. PreANSI(Pmatrix matrix)
  1669. PreANSI(Ppoint *point)
  1670. /*
  1671. ** \parambegin
  1672. ** \param{Pmatrix}{matrix}{3x3 matrix}{IN}
  1673. ** \param{Ppoint *}{point}{2D point}{IN}
  1674. ** \paramend
  1675. ** \blurb{This function transforms the 2D point \pardesc{point} by the
  1676. ** $3 \times 3$ matrix \pardesc{matrix}.}
  1677. */
  1678. {
  1679.   Ppoint temp;
  1680.   Pfloat w;
  1681.  
  1682.   temp.x = matrix[0][0] * point->x + matrix[0][1] * point->y + matrix[0][2];
  1683.   temp.y = matrix[1][0] * point->x + matrix[1][1] * point->y + matrix[1][2];
  1684.   return temp;
  1685. }  /* ptk_transform */
  1686.  
  1687. /*-------------------------------------------------------------------*/
  1688.  
  1689. /*function:external*/
  1690. extern void ptk_matrixtomatrix3(C(Pmatrix) mat, C(Pmatrix3) mat3)
  1691. PreANSI(Pmatrix mat)
  1692. PreANSI(Pmatrix mat3)
  1693. /*
  1694. ** \parambegin
  1695. ** \param{Pmatrix}{mat}{3x3 matrix}{IN}
  1696. ** \param{Pmatrix3}{mat3}{4x4 matrix}{OUT}
  1697. ** \paramend
  1698. ** \blurb{This function converts the $3 \times 3$ matrix \pardesc{mat}
  1699. ** to the $4 \times 4$ matrix \pardesc{mat3}, as follows:
  1700. ** $$  \left( \begin{array}{ccc}
  1701.            a & b & c\\
  1702.            d & e & f\\
  1703.            g & h & j
  1704.          \end{array} \right)
  1705.  
  1706. \rightarrow
  1707.  
  1708.  \left( \begin{array}{cccc}
  1709.            a & b & 0 & c\\
  1710.            d & e & 0 & f\\
  1711.        0 & 0 & 1 & 0\\
  1712.            g & h & 0 & j
  1713.          \end{array} \right) 
  1714. $$}
  1715. */
  1716. {
  1717.   ptk_unitmatrix3(mat3);
  1718.   mat3[0][0] = mat[0][0]; /* a */
  1719.   mat3[0][1] = mat[0][1]; /* b */
  1720.   mat3[1][0] = mat[1][0]; /* d */
  1721.   mat3[1][1] = mat[1][1]; /* e */
  1722.   mat3[0][3] = mat[0][2]; /* c */
  1723.   mat3[1][3] = mat[1][2]; /* f */
  1724.   mat3[3][0] = mat[2][0]; /* g */
  1725.   mat3[3][1] = mat[2][1]; /* h */
  1726.   mat3[3][3] = mat[2][2]; /* j */ 
  1727. }  /* ptk_matrixtomatrix3 */
  1728.  
  1729. /*----------------------------------------------------------------------------*/
  1730.  
  1731. /*function:external*/
  1732. extern void ptk_outputmatrix3(C(FILE *) fileptr, C(Pmatrix3) matrix, 
  1733.                               C(char *) string)
  1734. PreANSI(FILE *fileptr)
  1735. PreANSI(Pmatrix3 matrix)
  1736. PreANSI(char *string)
  1737. /*
  1738. ** \parambegin
  1739. ** \param{FILE *}{fileptr}{file pointer}{OUT}
  1740. ** \param{Pmatrix3}{matrix}{4x4 matrix}{IN}
  1741. ** \param{char *}{string}{character string}{IN}
  1742. ** \paramend
  1743. ** \blurb{This function outputs the $3\times 3$ matrix
  1744. **  \pardesc{matrix} and the message \pardesc{string}
  1745. ** to the file specified by \pardesc{fileptr}.}
  1746. */
  1747. {
  1748.   Pint ii;
  1749.  
  1750.   /* Output Title */
  1751.   fprintf(fileptr,"\n*************************************************\n");
  1752.   fprintf(fileptr,"%s\n", string);
  1753.   fprintf(fileptr,"           1          2          3          4\n");
  1754.   fprintf(fileptr,"-------------------------------------------------\n");
  1755.   for (ii = 0; ii <= 3; ii++)
  1756.     fprintf(fileptr,"%2ld | %11.4f%11.4f%11.4f%11.4f\n",
  1757.         ii + 1, matrix[ii][0], matrix[ii][1], matrix[ii][2], 
  1758.             matrix[ii][3]);
  1759.   fprintf(fileptr,"\n*************************************************\n");
  1760. }  /* ptk_outputmatrix3 */
  1761.  
  1762. /*--------------------------------------------------------------------------*/
  1763.  
  1764. static ptkboolean validbounds3(C(Plimit3 *) bound)
  1765. PreANSI(Plimit3 *bound)
  1766. /*
  1767. ** \parambegin
  1768. ** \param{Pint}{}{}{IN}
  1769. ** \paramend
  1770. ** description: checks if bounding box is sensibly defined (the min value
  1771. ** is less than max value in each axis).
  1772. ** input params: bound - 3D bounding box.
  1773. ** output params: none.
  1774. ** return value: TRUE if bounding box is valid, otherwise FALSE.
  1775. */
  1776. {
  1777.   return ((bound->x_min <= bound->x_max) && (bound->y_min <= bound->y_max) &&
  1778.       (bound->z_min <= bound->z_max));
  1779. }  /* validbounds3 */
  1780.  
  1781. /*--------------------------------------------------------------------------*/
  1782.  
  1783. static ptkboolean validbounds(C(Plimit *) bound)
  1784. PreANSI(Plimit *bound)
  1785. /*
  1786. ** \parambegin
  1787. ** \param{Pint}{}{}{IN}
  1788. ** \paramend
  1789. ** description: checks if bounding box is sensibly defined (the min value
  1790. ** is less than max value in each axis).
  1791. ** input params: bound - 3D bounding box.
  1792. ** output params: none.
  1793. ** return value: TRUE if bounding box is valid, otherwise FALSE.
  1794. */
  1795. {
  1796.   return ((bound->x_min <= bound->x_max) && (bound->y_min <= bound->y_max));
  1797. }  /* validbounds */
  1798.  
  1799. /*--------------------------------------------------------------------------*/
  1800.  
  1801. /*function:external*/
  1802. extern void ptk_box3tobox3(C(Plimit3 *) box1, C(Plimit3 *) box2, 
  1803.                            C(ptkboolean) preserve, C(Pcompose_type) operation, 
  1804.                            C(Pmatrix3) matrix, C(Pint *)error)
  1805. PreANSI(Plimit3 *box1)
  1806. PreANSI(Plimit3 *box2)
  1807. PreANSI(ptkboolean preserve)
  1808. PreANSI(Pcompose_type operation) 
  1809. PreANSI(Pmatrix3 matrix)
  1810. PreANSI(Pint *error)
  1811. /*
  1812. ** \parambegin
  1813. ** \param{Plimit3 *}{box1}{3D volume}{IN}
  1814. ** \param{Plimit3 *}{box2}{3D volume}{IN}
  1815. ** \param{ptkboolean}{preserve}{preserve aspect ratio}{IN}
  1816. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1817. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  1818. ** \param{Pint *}{error}{error code}{OUT}
  1819. ** \paramend
  1820. ** \blurb{This function computes a mapping from one 3D
  1821. ** box to another -- a 3D window to 3D viewport
  1822. ** transformation -- and concatenates this transformation
  1823. ** with \pardesc{matrix} on the
  1824. ** basis of \pardesc{operation}.If the parameters are invalid,
  1825. **  \pardesc{error} is set to 
  1826. ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
  1827. */
  1828. {
  1829.   Ppoint3 vpscale, winscale;
  1830.   Pmatrix3 temp;
  1831.   Ppoint3 winmin, winmax, vpmin, vpmax, wincentre, vpcentre;
  1832.   Pfloat minscale;
  1833.  
  1834.   *error = ptkcpcok;   /* assume OK */
  1835.   if (!validbounds3(box1)) 
  1836.   {
  1837.     *error = -1;
  1838.     return;
  1839.   }
  1840.   winmin = ptk_point3(box1->x_min, box1->y_min, box1->z_min);
  1841.   winmax = ptk_point3(box1->x_max, box1->y_max, box1->z_max);
  1842.   vpmin = ptk_point3(box2->x_min, box2->y_min, box2->z_min);
  1843.   vpmax = ptk_point3(box2->x_max, box2->y_max, box2->z_max);
  1844.  
  1845.   winscale = ptk_subv3(&winmax, &winmin);
  1846.   wincentre = ptk_scalev3(&winscale, 0.5);
  1847.   wincentre = ptk_addv3(&winmin, &wincentre);
  1848.   wincentre = ptk_scalev3(&wincentre, -1.0);
  1849.   vpscale = ptk_subv3(&vpmax, &vpmin);
  1850.   vpcentre = ptk_scalev3(&vpscale, 0.5);
  1851.   vpcentre = ptk_addv3(&vpmin, &vpcentre);
  1852.   if (winscale.x == 0.0)
  1853.     winscale.x = 1.0;
  1854.   else
  1855.     winscale.x = vpscale.x / winscale.x;
  1856.   if (winscale.y == 0.0)
  1857.     winscale.y = 1.0;
  1858.   else
  1859.     winscale.y = vpscale.y / winscale.y;
  1860.   if (winscale.z == 0.0)
  1861.     winscale.z = 1.0;
  1862.   else
  1863.     winscale.z = vpscale.z / winscale.z;
  1864.   if (preserve)
  1865.   {
  1866.     minscale = PTKMIN(winscale.x, winscale.y);
  1867.     minscale = PTKMIN(minscale, winscale.z);
  1868.     winscale = ptk_point3(minscale, minscale, minscale);
  1869.   }
  1870.   ptk_shift3(&vpcentre, PTYPE_REPLACE, temp);
  1871.   ptk_scale3(&winscale, PTYPE_PRECONCAT, temp);
  1872.   ptk_shift3(&wincentre, PTYPE_PRECONCAT, temp);
  1873.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  1874. }  /* ptk_box3tobox3 */
  1875.  
  1876. /*--------------------------------------------------------------------------*/
  1877.  
  1878. /*function:external*/
  1879. extern void ptk_boxtobox(C(Plimit *) box1, C(Plimit *) box2, 
  1880.                            C(ptkboolean) preserve, C(Pcompose_type) operation, 
  1881.                            C(Pmatrix) matrix, C(Pint *)error)
  1882. PreANSI(Plimit *box1)
  1883. PreANSI(Plimit *box2)
  1884. PreANSI(ptkboolean preserve)
  1885. PreANSI(Pcompose_type operation) 
  1886. PreANSI(Pmatrix matrix)
  1887. PreANSI(Pint *error)
  1888. /*
  1889. ** \parambegin
  1890. ** \param{Plimit *}{box1}{2D box}{IN}
  1891. ** \param{Plimit *}{box2}{2D box}{IN}
  1892. ** \param{ptkboolean}{preserve}{preserve aspect ratio}{IN}
  1893. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1894. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  1895. ** \param{Pint *}{error}{error code}{OUT}
  1896. ** \paramend
  1897. ** \blurb{This function computes
  1898. **  a mapping from one 2D box to another -- a 2D window to 2D viewport
  1899. ** transformation ---and concatenates this transformation
  1900. ** it with \pardesc{matrix} on the
  1901. ** basis of \pardesc{operation}.
  1902. ** If the parameters are invalid, \pardesc{error} is set to 
  1903. ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
  1904. */
  1905. {
  1906.   Ppoint vpscale, winscale;
  1907.   Pmatrix temp;
  1908.   Ppoint winmin, winmax, vpmin, vpmax, wincentre, vpcentre;
  1909.   Pfloat minscale;
  1910.  
  1911.   *error = ptkcpcok;   /* assume OK */
  1912.   if (!validbounds(box1)) 
  1913.   {
  1914.     *error = -1;
  1915.     return;
  1916.   }
  1917.   winmin = ptk_point(box1->x_min, box1->y_min);
  1918.   winmax = ptk_point(box1->x_max, box1->y_max);
  1919.   vpmin = ptk_point(box2->x_min, box2->y_min);
  1920.   vpmax = ptk_point(box2->x_max, box2->y_max);
  1921.  
  1922.   winscale = ptk_subv(&winmax, &winmin);
  1923.   wincentre = ptk_scalev(&winscale, 0.5);
  1924.   wincentre = ptk_addv(&winmin, &wincentre);
  1925.   wincentre = ptk_scalev(&wincentre, -1.0);
  1926.   vpscale = ptk_subv(&vpmax, &vpmin);
  1927.   vpcentre = ptk_scalev(&vpscale, 0.5);
  1928.   vpcentre = ptk_addv(&vpmin, &vpcentre);
  1929.   if (winscale.x == 0.0)
  1930.     winscale.x = 1.0;
  1931.   else
  1932.     winscale.x = vpscale.x / winscale.x;
  1933.   if (winscale.y == 0.0)
  1934.     winscale.y = 1.0;
  1935.   else
  1936.     winscale.y = vpscale.y / winscale.y;
  1937.   if (preserve)
  1938.   {
  1939.     minscale = PTKMIN(winscale.x, winscale.y);
  1940.     winscale = ptk_point(minscale, minscale);
  1941.   }
  1942.   ptk_shift(&vpcentre, PTYPE_REPLACE, temp);
  1943.   ptk_scale(&winscale, PTYPE_PRECONCAT, temp);
  1944.   ptk_shift(&wincentre, PTYPE_PRECONCAT, temp);
  1945.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  1946. }  /* ptk_boxtobox */
  1947.  
  1948. /*--------------------------------------------------------------------------*/
  1949.  
  1950. /*function:external*/
  1951. extern void ptk_accumulatetran3(C(Ppoint3 *) fixed, C(Ppoint3 *) shift, 
  1952.                                 C(Pfloat) rotx, C(Pfloat) roty, 
  1953.                                 C(Pfloat) rotz, C(Ppoint3 *) scale, 
  1954.                                 C(Pcompose_type) operation, 
  1955.                                 C(Pmatrix3) matrix)
  1956. PreANSI(Ppoint3 *fixed)
  1957. PreANSI(Ppoint3 *shift) 
  1958. PreANSI(Pfloat rotx)
  1959. PreANSI(Pfloat roty)
  1960. PreANSI(Pfloat rotz) 
  1961. PreANSI(Ppoint3 *scale)
  1962. PreANSI(Pcompose_type operation) 
  1963. PreANSI(Pmatrix3 matrix)
  1964. /*
  1965. ** \parambegin
  1966. ** \param{Ppoint3 *}{fixed}{origin}{IN}
  1967. ** \param{Ppoint3 *}{shift}{shift factor}{IN}
  1968. ** \param{Pfloat}{rotx}{x rotation}{IN}
  1969. ** \param{Pfloat}{rotx}{y rotation}{IN}
  1970. ** \param{Pfloat}{rotx}{z rotation}{IN}
  1971. ** \param{Ppoint3 *}{scale}{scale factor}{IN}
  1972. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  1973. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  1974. ** \paramend
  1975. ** \blurb{This function computes the specified 3D
  1976. ** shift, scale and rotate transformations, in the order 
  1977. ** scale, rotate, shift, 
  1978. ** and then concatenates the resulting transformation
  1979. ** with the specified matrix on the basis of \pardesc{operation}.}
  1980. */
  1981. {
  1982.   Pmatrix3 temp;
  1983.   Ppoint3 mfixed;
  1984.  
  1985.   mfixed = ptk_scalev3(fixed, -1.0);
  1986.   ptk_shift3(shift, PTYPE_REPLACE, temp);
  1987.   ptk_shift3(fixed, PTYPE_PRECONCAT, temp);
  1988.   ptk_rotate3(rotx, PTKEXAXIS, PTYPE_PRECONCAT, temp);
  1989.   ptk_rotate3(roty, PTKEYAXIS, PTYPE_PRECONCAT, temp);
  1990.   ptk_rotate3(rotz, PTKEZAXIS, PTYPE_PRECONCAT, temp);
  1991.   ptk_scale3(scale, PTYPE_PRECONCAT, temp);
  1992.   ptk_shift3(&mfixed, PTYPE_PRECONCAT, temp);
  1993.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  1994. }  /* ptk_accumulatetran3 */
  1995.  
  1996. /*--------------------------------------------------------------------------*/
  1997.  
  1998. /*function:external*/
  1999. extern void ptk_accumulatetran(C(Ppoint *) fixed, C(Ppoint *) shift, 
  2000.                                 C(Pfloat) rot, C(Ppoint *) scale, 
  2001.                                 C(Pcompose_type) operation, 
  2002.                                 C(Pmatrix) matrix)
  2003. PreANSI(Ppoint *fixed)
  2004. PreANSI(Ppoint *shift) 
  2005. PreANSI(Pfloat rot)
  2006. PreANSI(Ppoint *scale)
  2007. PreANSI(Pcompose_type operation) 
  2008. PreANSI(Pmatrix matrix)
  2009. /*
  2010. ** \parambegin
  2011. ** \param{Ppoint *}{fixed}{origin}{IN}
  2012. ** \param{Ppoint *}{shift}{shift factor}{IN}
  2013. ** \param{Pfloat}{rotx}{x rotation}{IN}
  2014. ** \param{Ppoint *}{scale}{scale factor}{IN}
  2015. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  2016. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  2017. ** \paramend
  2018. ** \blurb{This function computes the specified 2D
  2019. ** shift, scale and rotate transformations, in the order 
  2020. ** scale, rotate, shift, 
  2021. ** and then concatenates the resulting transformation
  2022. ** with the specified matrix on the basis of \pardesc{operation}.}
  2023. */
  2024. {
  2025.   Pmatrix temp;
  2026.   Ppoint mfixed;
  2027.  
  2028.   mfixed = ptk_scalev(fixed, -1.0);
  2029.   ptk_shift(shift, PTYPE_REPLACE, temp);
  2030.   ptk_shift(fixed, PTYPE_PRECONCAT, temp);
  2031.   ptk_rotate(rot, PTKEXAXIS, PTYPE_PRECONCAT, temp);
  2032.   ptk_scale(scale, PTYPE_PRECONCAT, temp);
  2033.   ptk_shift(&mfixed, PTYPE_PRECONCAT, temp);
  2034.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  2035. }  /* ptk_accumulatetran */
  2036.  
  2037. /*--------------------------------------------------------------------------*/
  2038.  
  2039. /*function:external*/
  2040. extern void ptk_evalvieworientation3(C(Ppoint3 *) viewrefpoint, 
  2041.                                     C(Ppoint3 *) viewplanenormal,
  2042.                     C(Ppoint3 *) viewupvector, 
  2043.                                     C(Pcompose_type) operation, 
  2044.                                     C(Pmatrix3) matrix, C(Pint *)error)
  2045. PreANSI(Ppoint3 *viewrefpoint)
  2046. PreANSI(Ppoint3 *viewplanenormal)
  2047. PreANSI(Ppoint3 *viewupvector)
  2048. PreANSI(Pcompose_type operation) 
  2049. PreANSI(Pmatrix3 matrix)
  2050. PreANSI(Pint *error)
  2051. /*
  2052. ** \parambegin
  2053. ** \param{Ppoint3 *}{viewrefpoint}{view reference point}{IN}
  2054. ** \param{Ppoint3 *}{viewplanenormal}{view plane normal}{IN}
  2055. ** \param{Ppoint3 *}{viewupvector}{view up vector}{IN}
  2056. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  2057. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  2058. ** \param{Pint *}{error}{error code}{OUT}
  2059. ** \paramend
  2060. ** \blurb{This function computes 
  2061. ** a 3D PHIGS view orientation matrix on the basis of
  2062. ** a specified view reference point (\pardesc{viewrefpoint}), a
  2063. ** view plane normal (\pardesc{viewplanenormal}) and a view up vector
  2064. ** (\pardesc{viewupvector}). If the function succeeds,
  2065. **  \pardesc{error} is set to 
  2066. **  \pardesc{ptkcpcok}. Otherwise,
  2067. ** \pardesc{error} is 61 if the view plane normal is null,
  2068. ** 63 if the view up vector is null,
  2069. ** and 58 if the cross product of the view up vector
  2070. ** and the view plane normal is null.}
  2071. */
  2072. {
  2073.   Ppoint3 uaxis, vaxis, naxis, vuvxvpn;
  2074.   Pmatrix3 temp;
  2075.  
  2076.   *error = ptkcpcok;
  2077.   if (ptk_nullv3(viewplanenormal)) 
  2078.   {
  2079.     *error = 61;
  2080.   } 
  2081.   else 
  2082.   if (ptk_nullv3(viewupvector)) 
  2083.   {
  2084.     *error = 63;
  2085.   } 
  2086.   else 
  2087.   {
  2088.     vuvxvpn = ptk_crossv3(viewupvector, viewplanenormal);
  2089.     if (ptk_nullv3(&vuvxvpn))
  2090.       *error = 58;
  2091.     else
  2092.       *error = ptkcpcok;
  2093.   }
  2094.   if (*error != ptkcpcok) 
  2095.   {
  2096.     return;
  2097.   }  /* if error = 0 */
  2098.   /* generate from viewplanenormal and viewupvector a set
  2099.   ** of orthogonal axes denoted by uaxis, vaxis, naxis 
  2100.   */
  2101.   naxis = ptk_unitv3(viewplanenormal);
  2102.   uaxis = ptk_unitv3(&vuvxvpn);
  2103.   vaxis = ptk_crossv3(&naxis, &uaxis);
  2104.   vaxis = ptk_unitv3(&vaxis);
  2105.   /* form view orientation matrix */
  2106.   ptk_unitmatrix3(temp);
  2107.   temp[0][0] = uaxis.x;
  2108.   temp[0][1] = uaxis.y;
  2109.   temp[0][2] = uaxis.z;
  2110.   temp[0][3] = -(uaxis.x * viewrefpoint->x + uaxis.y * 
  2111.                    viewrefpoint->y + uaxis.z * viewrefpoint->z);
  2112.   temp[1][0] = vaxis.x;
  2113.   temp[1][1] = vaxis.y;
  2114.   temp[1][2] = vaxis.z;
  2115.   temp[1][3] = -(vaxis.x * viewrefpoint->x + vaxis.y * 
  2116.                    viewrefpoint->y + vaxis.z * viewrefpoint->z);
  2117.   temp[2][0] = naxis.x;
  2118.   temp[2][1] = naxis.y;
  2119.   temp[2][2] = naxis.z;
  2120.   temp[2][3] = -(naxis.x * viewrefpoint->x + naxis.y * 
  2121.                    viewrefpoint->y + naxis.z * viewrefpoint->z);
  2122.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  2123. }  /* ptk_evalvieworientation3 */
  2124.  
  2125. /*--------------------------------------------------------------------------*/
  2126.  
  2127. /*function:external*/
  2128. extern void ptk_evalvieworientation(C(Ppoint *) viewrefpoint, 
  2129.                     C(Ppoint *) viewupvector, 
  2130.                                     C(Pcompose_type) operation, 
  2131.                                     C(Pmatrix) matrix, C(Pint *)error)
  2132. PreANSI(Ppoint *viewrefpoint)
  2133. PreANSI(Ppoint *viewupvector)
  2134. PreANSI(Pcompose_type operation) 
  2135. PreANSI(Pmatrix matrix)
  2136. PreANSI(Pint *error)
  2137. /*
  2138. ** \parambegin
  2139. ** \param{Ppoint *}{viewrefpoint}{view reference point}{IN}
  2140. ** \param{Ppoint *}{viewupvector}{view up vector}{IN}
  2141. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  2142. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  2143. ** \param{Pint *}{error}{error code}{OUT}
  2144. ** \paramend
  2145. ** \blurb{This function computes 
  2146. ** a 2D PHIGS view orientation matrix on the basis of
  2147. ** a specified view reference point (\pardesc{viewrefpoint})
  2148. ** and a view up vector
  2149. ** (\pardesc{viewupvector}).  If the function succeeds,
  2150. ** \pardesc{error} is set to \pardesc{ptkcpcok}. Otherwise, 
  2151. ** \pardesc{error} is 63 if the view up vector is null.}
  2152. */
  2153. {
  2154.   Ppoint3 uaxis, vaxis, naxis, vuvxvpn;
  2155.   Pmatrix temp;
  2156.  
  2157.   *error = ptkcpcok;
  2158.   if (ptk_nullv(viewupvector)) 
  2159.   {
  2160.     *error = 63;
  2161.   } 
  2162.   if (*error != ptkcpcok) 
  2163.   {
  2164.     return;
  2165.   }  /* if error = 0 */
  2166.   /* generate from viewplanenormal and viewupvector a set
  2167.   ** of orthogonal axes denoted by uaxis, vaxis, naxis 
  2168.   */
  2169.   naxis = ptk_point3(0.0, 0.0, 1.0);
  2170.   vaxis = ptk_point3(viewupvector->x, viewupvector->y, 0.0);
  2171.   vaxis = ptk_unitv3(&vaxis);
  2172. /*  uaxis = ptk_crossv3(naxis, vaxis); */
  2173.   uaxis = ptk_unitv3(&uaxis);
  2174.   /* form view orientation matrix */
  2175.   ptk_unitmatrix(temp);
  2176.   temp[0][0] = uaxis.x;
  2177.   temp[0][1] = uaxis.y;
  2178.   temp[0][2] = uaxis.z;
  2179. /*  temp[0][3] = -(uaxis.x * viewrefpoint->x + uaxis.y * viewrefpoint->y); */
  2180.   temp[1][0] = vaxis.x;
  2181.   temp[1][1] = vaxis.y;
  2182.   temp[1][2] = 0.0;
  2183. /*  temp[1][3] = -(vaxis.x * viewrefpoint->x + vaxis.y * viewrefpoint->y); */
  2184.   temp[2][0] = naxis.x;
  2185.   temp[2][1] = naxis.y;
  2186.   temp[2][2] = naxis.z;
  2187.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  2188. }  /* ptk_evalvieworientation */
  2189.  
  2190. /*--------------------------------------------------------------------------*/
  2191.  
  2192. /*function:external*/
  2193. extern void ptk_evalviewmapping3(C(Plimit3 *) wlimits, C(Plimit3 *) vlimits, 
  2194.                                 C(Pproj_type) viewtype, 
  2195.                                 C(Ppoint3 *) refpoint, C(Pfloat) vplanedist, 
  2196.                                 C(Pcompose_type) operation, C(Pmatrix3) matrix, 
  2197.                                 C(Pint *)error)
  2198. PreANSI(Plimit3 *wlimits)
  2199. PreANSI(Plimit3 *vlimits) 
  2200. PreANSI(Pproj_type viewtype) 
  2201. PreANSI(Ppoint3 *refpoint)
  2202. PreANSI(Pfloat vplanedist) 
  2203. PreANSI(Pcompose_type operation)
  2204. PreANSI(Pmatrix3 matrix) 
  2205. PreANSI(Pint *error)
  2206. /*
  2207. ** \parambegin
  2208. ** \param{Plimit3 *}{wlimits}{window limits}{IN}
  2209. ** \param{Plimit3 *}{vlimits}{viewport limits}{IN}
  2210. ** \param{Pproj\_type}{viewtype}{projection type}{IN}
  2211. ** \param{Ppoint3 *}{refpoint}{projection reference point}{IN}
  2212. ** \param{Pfloat}{vplanedist}{view plane distance}{IN}
  2213. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  2214. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  2215. ** \param{Pint *}{error}{error code}{OUT}
  2216. ** \paramend
  2217. ** \blurb{This function evaluates a 3D PHIGS view mapping matrix.
  2218. ** If the function succeeds,
  2219. ** \pardesc{error} is set to \pardesc{ptkcpcok}. Otherwise, 
  2220. ** \pardesc{error} is 
  2221. ** 329 if the window limits are not valid,
  2222. ** 336 if the back plane is in front of front plane,
  2223. ** 330 if the viewport limits are not valid,
  2224. ** 335 if the projection reference point is on the view plane,
  2225. ** and 340 if the projection reference point is between front and back planes.}
  2226. */
  2227. {
  2228.   Ppoint3 projectordirn, viewwindowcentre, vref;
  2229.   Pmatrix3 temp;
  2230.   Plimit3 templim;
  2231.   Pint err;
  2232.  
  2233.   templim = *wlimits;
  2234.   *error = ptkcpcok;
  2235.   /* check for valid window limits */
  2236.   if (templim.x_min >= templim.x_max ||
  2237.       templim.y_min >= templim.y_max) 
  2238.   {
  2239.     *error = 329;
  2240.     return;
  2241.   }
  2242.   /* check that back plane is not in front of front plane */
  2243.   if (templim.z_min >= templim.z_max) 
  2244.   {
  2245.     *error = 336;
  2246.     return;
  2247.   }
  2248.   /* check for valid viewport limits */
  2249.   if (!validbounds3(vlimits)) 
  2250.   {
  2251.     *error = 330;
  2252.     return;
  2253.   }
  2254.   /* check that PRP not positioned on the view plane */
  2255.   if (ptk_equal(refpoint->z, vplanedist)) 
  2256.   {
  2257.     *error = 335;
  2258.     return;
  2259.   }
  2260.   /* check that PRP not between the front and back planes */
  2261.   if (refpoint->z < templim.z_max && refpoint->z > templim.z_min) 
  2262.   {
  2263.     *error = 340;
  2264.     return;
  2265.   }
  2266.   /* if everything is OK error is set to zero and a value computed for 
  2267.   ** matrix.
  2268.   */
  2269.   *error = ptkcpcok;
  2270.   ptk_unitmatrix3(temp);
  2271.   switch (viewtype) 
  2272.   {
  2273.   case PTYPE_PERSPECT:
  2274.     temp[0][0] = refpoint->z - vplanedist;
  2275.     temp[0][2] = -refpoint->x;
  2276.     temp[0][3] = refpoint->x * vplanedist;
  2277.     temp[1][0] = 0.0;
  2278.     temp[1][1] = refpoint->z - vplanedist;
  2279.     temp[1][2] = -refpoint->y;
  2280.     temp[1][3] = refpoint->y * vplanedist;
  2281.     temp[2][2] = 0.0;
  2282.     temp[2][3] = 1.0;
  2283.     temp[3][2] = -1.0;
  2284.     temp[3][3] = refpoint->z;
  2285.     break;
  2286.  
  2287.   case PTYPE_PARAL:
  2288.     viewwindowcentre.x = (templim.x_min + templim.x_max) / 2;
  2289.     viewwindowcentre.y = (templim.y_min + templim.y_max) / 2;
  2290.     viewwindowcentre.z = vplanedist;
  2291.     projectordirn = ptk_subv3(refpoint, &viewwindowcentre);
  2292.     projectordirn = ptk_scalev3(&projectordirn, 
  2293.                                1.0 / (refpoint->z - vplanedist));
  2294.     temp[0][2] = -projectordirn.x;
  2295.     temp[0][3] = vplanedist * projectordirn.x;
  2296.     temp[1][2] = -projectordirn.y;
  2297.     temp[1][3] = vplanedist * projectordirn.y;
  2298.     break;
  2299.   }
  2300.   /* The z coordinate is unchanged by the transformation matrix above for
  2301.   ** parallel projection types and changed for perspective projection
  2302.   ** types. Different action is therefore required in each case so that
  2303.   ** points at the front and back planes are always mapped onto points at
  2304.   ** the front and back of the projection viewport 
  2305.   */
  2306.   if (viewtype == PTYPE_PERSPECT) 
  2307.   {
  2308.     templim.z_min = 1.0 / (refpoint->z - templim.z_min);
  2309.     templim.z_max = 1.0 / (refpoint->z - templim.z_max);
  2310.   }
  2311.   ptk_box3tobox3(&templim, vlimits, FALSE, PTYPE_POSTCONCAT, temp, &err);
  2312.   ptk_concatenatematrix3(operation, temp, matrix, matrix);
  2313. }  /* ptk_evalviewmapping3 */
  2314.  
  2315. /*--------------------------------------------------------------------------*/
  2316.  
  2317. /*function:external*/
  2318. extern void ptk_evalviewmapping(C(Plimit *) wlimits, C(Plimit *) vlimits, 
  2319.                                 C(Pcompose_type) operation, C(Pmatrix) matrix, 
  2320.                                 C(Pint *)error)
  2321. PreANSI(Plimit *wlimits)
  2322. PreANSI(Plimit *vlimits) 
  2323. PreANSI(Pcompose_type operation)
  2324. PreANSI(Pmatrix matrix) 
  2325. PreANSI(Pint *error)
  2326. /*
  2327. ** \parambegin
  2328. ** \param{Plimit *}{wlimits}{window limits}{IN}
  2329. ** \param{Plimit *}{vlimits}{viewport limits}{IN}
  2330. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  2331. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  2332. ** \param{Pint *}{error}{error code}{OUT}
  2333. ** \paramend
  2334. ** \blurb{This function evaluates a 2d PHIGS view mapping matrix.
  2335. ** If the function succeeds,
  2336. ** \pardesc{error} is set to \pardesc{ptkcpcok}. Otherwise, 
  2337. ** \pardesc{error} is
  2338. ** 329 if the window limits are not valid,
  2339. ** and 330 if the viewport limits are not valid.}
  2340. */
  2341. {
  2342.   Pmatrix temp;
  2343.   Pint err;
  2344.  
  2345.   *error = ptkcpcok;
  2346.   /* check for valid window limits */
  2347.   if (wlimits->x_min >= wlimits->x_max ||
  2348.       wlimits->y_min >= wlimits->y_max) 
  2349.   {
  2350.     *error = 329;
  2351.     return;
  2352.   }
  2353.   /* check for valid viewport limits */
  2354.   if (!validbounds(vlimits)) 
  2355.   {
  2356.     *error = 330;
  2357.     return;
  2358.   }
  2359.   /* if everything is OK error is set to zero and a value computed for 
  2360.   ** matrix.
  2361.   */
  2362.   *error = ptkcpcok;
  2363.   ptk_unitmatrix(temp);
  2364.   ptk_boxtobox(wlimits, vlimits, FALSE, PTYPE_POSTCONCAT, temp, &err);
  2365.   ptk_concatenatematrix(operation, temp, matrix, matrix);
  2366. }  /* ptk_evalviewmapping */
  2367.  
  2368. /*--------------------------------------------------------------------------*/
  2369.  
  2370. /*function:external*/
  2371. extern void ptk_stackmatrix3(C(Pmatrix3) matrix)
  2372. PreANSI(Pmatrix3 matrix)
  2373. /*
  2374. ** \parambegin
  2375. ** \param{Pmatrix3}{matrix}{4x4 matrix}{IN}
  2376. ** \paramend
  2377. ** \blurb{This function pushes
  2378. ** the $4 \times 4$ matrix \pardesc{matrix} onto the transformation stack.}
  2379. */
  2380. {
  2381.   newone3 = (ptksstack3 *)malloc(sizeof(ptksstack3));
  2382.   newone3->ptkpnext = listhead3;
  2383.   memcpy(newone3->ptkamat, matrix, sizeof(Pmatrix3));
  2384.   listhead3 = newone3;
  2385. }  /* ptk_stackmatrix3 */
  2386.  
  2387. /*--------------------------------------------------------------------------*/
  2388.  
  2389. /*function:external*/
  2390. extern void ptk_stackmatrix(C(Pmatrix) matrix)
  2391. PreANSI(Pmatrix matrix)
  2392. /*
  2393. ** \parambegin
  2394. ** \param{Pmatrix}{matrix}{3x3 matrix}{IN}
  2395. ** \paramend
  2396. ** \blurb{This function pushes the $3 \times 3$ matrix \pardesc{matrix}
  2397. **  onto the transformation stack.}
  2398. */
  2399. {
  2400.   newone = (ptksstack *)malloc(sizeof(ptksstack));
  2401.   newone->ptkpnext = listhead;
  2402.   memcpy(newone->ptkamat, matrix, sizeof(Pmatrix));
  2403.   listhead = newone;
  2404. }  /* ptk_stackmatrix */
  2405.  
  2406. /*--------------------------------------------------------------------------*/
  2407.  
  2408. /*function:external*/
  2409. extern void ptk_unstackmatrix3(C(Pmatrix3) matrix)
  2410. PreANSI(Pmatrix3 matrix)
  2411. /*
  2412. ** \parambegin
  2413. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  2414. ** \paramend
  2415. ** \blurb{This function
  2416. ** pops a $4 \times 4$  matrix 
  2417. ** from the transformation stack and returns it in
  2418. ** \pardesc{matrix}.}
  2419. */
  2420. {
  2421.   if (listhead3 == NULL)
  2422.     return;
  2423.   memcpy(matrix, listhead3->ptkamat, sizeof(Pmatrix3));
  2424.   newone3 = listhead3->ptkpnext;
  2425.   free(listhead3);
  2426.   listhead3 = newone3;
  2427. }  /* ptk_unstackmatrix3 */
  2428.  
  2429. /*--------------------------------------------------------------------------*/
  2430.  
  2431. /*function:external*/
  2432. extern void ptk_unstackmatrix(C(Pmatrix) matrix)
  2433. PreANSI(Pmatrix matrix)
  2434. /*
  2435. ** \parambegin
  2436. ** \param{Pmatrix}{matrix}{3x3 matrix}{OUT}
  2437. ** \paramend
  2438. ** \blurb{This function pops a $3 \times 3$ matrix
  2439. ** from the transformation stack and returns it in
  2440. ** \pardesc{matrix}.}
  2441. */
  2442. {
  2443.   if (listhead == NULL)
  2444.     return;
  2445.   memcpy(matrix, listhead->ptkamat, sizeof(Pmatrix));
  2446.   newone = listhead->ptkpnext;
  2447.   free(listhead);
  2448.   listhead = newone;
  2449. }  /* ptk_unstackmatrix */
  2450.  
  2451. /*--------------------------------------------------------------------------*/
  2452.  
  2453. /*function:external*/
  2454. extern void ptk_examinestackmatrix3(C(Pmatrix3) matrix)
  2455. PreANSI(Pmatrix3 matrix)
  2456. /*
  2457. ** \parambegin
  2458. ** \param{Pmatrix3}{matrix}{4x4 matrix}{IN}
  2459. ** \paramend
  2460. ** \blurb{This function returns the top entry on the transformation stack.
  2461. ** The stack is not disturbed.}
  2462. */
  2463. {
  2464.   if (listhead3 != NULL)
  2465.     memcpy(matrix, listhead3->ptkamat, sizeof(Pmatrix3));
  2466. }  /* ptk_examinestackmatrix3 */
  2467.  
  2468. /*--------------------------------------------------------------------------*/
  2469.  
  2470. /*function:external*/
  2471. extern void ptk_examinestackmatrix(C(Pmatrix) matrix)
  2472. PreANSI(Pmatrix matrix)
  2473. /*
  2474. ** \parambegin
  2475. ** \param{Pmatrix}{matrix}{3x3 matrix}{IN}
  2476. ** \paramend
  2477. ** \blurb{This function  returns the top entry on the transformation stack.
  2478. ** The stack is not disturbed.}
  2479. */
  2480. {
  2481.   if (listhead != NULL)
  2482.     memcpy(matrix, listhead->ptkamat, sizeof(Pmatrix));
  2483. }  /* ptk_examinestackmatrix */
  2484.  
  2485. /*--------------------------------------------------------------------------*/
  2486.  
  2487. /*function:external*/
  2488. extern void ptk_3ptto3pt(C(Ppoint3 *) p1, C(Ppoint3 *) p2, C(Ppoint3 *) p3, 
  2489.                          C(Ppoint3 *) q1, C(Ppoint3 *) q2, C(Ppoint3 *) q3, 
  2490.                          C(Pcompose_type) operation, C(Pmatrix3) matrix, 
  2491.                          C(Pint *)error)
  2492. PreANSI(Ppoint3 *p1)
  2493. PreANSI(Ppoint3 *p2)
  2494. PreANSI(Ppoint3 *p3)
  2495. PreANSI(Ppoint3 *q1) 
  2496. PreANSI(Ppoint3 *q2)
  2497. PreANSI(Ppoint3 *q3)
  2498. PreANSI(Pcompose_type operation) 
  2499. PreANSI(Pmatrix3 matrix)
  2500. PreANSI(Pint *error)
  2501. /*
  2502. ** \parambegin
  2503. ** \param{Ppoint3 *}{p1}{3D point}{IN}
  2504. ** \param{Ppoint3 *}{p2}{3D point}{IN}
  2505. ** \param{Ppoint3 *}{p3}{3D point}{IN}
  2506. ** \param{Ppoint3 *}{q1}{3D point}{IN}
  2507. ** \param{Ppoint3 *}{q2}{3D point}{IN}
  2508. ** \param{Ppoint3 *}{q3}{3D point}{IN}
  2509. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  2510. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  2511. ** \param{Pint *}{error}{error code}{OUT}
  2512. ** \paramend
  2513. ** \blurb{This function returns the 3 point to 3 point transformation as
  2514. ** described in \cite{mort:geom}, pages 353--355.
  2515. ** The transformation has the following properties:
  2516. ** \pardesc{p1} is transformed onto \pardesc{q1};
  2517. ** the vector \pardesc{(p2-p1)} is transformed to be parallel to the vector 
  2518. ** \pardesc{(q2-q1)};
  2519. ** the plane containing the three points \pardesc{p1, p2, p3} is
  2520. ** transformed into  the plane containing \pardesc{q1, q2, q3}.
  2521. ** The transformation is concatenated with the $4 \times 4$ matrix
  2522. ** \pardesc{matrix} on the basis of \pardesc{operation}.
  2523. ** If the parameters are invalid, \pardesc{error} is set to 
  2524. ** -1. Otherwise, its value is \pardesc{ptkcpcok}.}
  2525. */
  2526. {
  2527.   Pmatrix3 tempmat, invvmatrix, wmatrix;
  2528.   Ppoint3 v1, v2, v3, w1, w2, w3, temp;
  2529.  
  2530.   /* form 2 right-hand orthogonal systems from the given points */
  2531.   *error = ptkcpcok;
  2532.   v1 = ptk_subv3(p2, p1);
  2533.   v1 = ptk_unitv3(&v1);
  2534.   v3 = ptk_subv3(p3, p1);
  2535.   v3 = ptk_crossv3(&v1, &v3);
  2536.   v3 = ptk_unitv3(&v3);
  2537.   v2 = ptk_crossv3(&v3, &v1);
  2538.   v2 = ptk_unitv3(&v2);
  2539.   temp = ptk_subv3(q2, q1);
  2540.   w1 = ptk_unitv3(&temp);
  2541.   temp = ptk_subv3(q3, q1);
  2542.   temp = ptk_crossv3(&w1, &temp);
  2543.   w3 = ptk_unitv3(&temp);
  2544.   temp = ptk_crossv3(&w3, &w1);
  2545.   w2 = ptk_unitv3(&temp);
  2546.  
  2547.   /* if any of v1..v3, w1..w3 are null vectors then the parameters are 
  2548.   ** invalid.
  2549.   */
  2550.   if (ptk_nullv3(&v1) || ptk_nullv3(&v2) || ptk_nullv3(&v3) || 
  2551.       ptk_nullv3(&w1) || ptk_nullv3(&w2) || ptk_nullv3(&w3)) 
  2552.   {
  2553.     *error = -1;
  2554.     return;
  2555.   }
  2556.   *error = ptkcpcok;
  2557.   ptk_unitmatrix3(invvmatrix);
  2558.   invvmatrix[0][0] = v1.x;
  2559.   invvmatrix[0][1] = v1.y;
  2560.   invvmatrix[0][2] = v1.z;
  2561.   invvmatrix[1][0] = v2.x;
  2562.   invvmatrix[1][1] = v2.y;
  2563.   invvmatrix[1][2] = v2.z;
  2564.   invvmatrix[2][0] = v3.x;
  2565.   invvmatrix[2][1] = v3.y;
  2566.   invvmatrix[2][2] = v3.z;
  2567.   ptk_unitmatrix3(wmatrix);
  2568.   wmatrix[0][0] = w1.x;
  2569.   wmatrix[1][0] = w1.y;
  2570.   wmatrix[2][0] = w1.z;
  2571.   wmatrix[0][1] = w2.x;
  2572.   wmatrix[1][1] = w2.y;
  2573.   wmatrix[2][1] = w2.z;
  2574.   wmatrix[0][2] = w3.x;
  2575.   wmatrix[1][2] = w3.y;
  2576.   wmatrix[2][2] = w3.z;
  2577.  
  2578.   /* shift points at origin onto q1 */
  2579.  
  2580.   ptk_shift3(q1, PTYPE_REPLACE, tempmat);
  2581.   ptk_concatenatematrix3(PTYPE_PRECONCAT, wmatrix, wmatrix, tempmat);
  2582.   ptk_concatenatematrix3(PTYPE_PRECONCAT, invvmatrix, tempmat, tempmat);
  2583.  
  2584.   /* shift p1 to origin */
  2585.   temp = ptk_scalev3(p1, -1.0);
  2586.   ptk_shift3(&temp, PTYPE_PRECONCAT, tempmat);
  2587.   ptk_concatenatematrix3(operation, tempmat, matrix, matrix);
  2588. }  /* ptk_3ptto3pt */
  2589.  
  2590. /*--------------------------------------------------------------------------*/
  2591.  
  2592. /*function:external*/
  2593. extern void ptk_0to3pt(C(Ppoint3 *) origin, C(Ppoint3 *) xdirn, 
  2594.               C(Ppoint3 *) ydirn, C(Pcompose_type) operation, C(Pmatrix3) matrix)
  2595. PreANSI(Ppoint3 *origin)
  2596. PreANSI(Ppoint3 *xdirn)
  2597. PreANSI(Ppoint3 *ydirn) 
  2598. PreANSI(Pcompose_type operation)
  2599. PreANSI(Pmatrix3 matrix)
  2600. /*
  2601. ** \parambegin
  2602. ** \param{Ppoint3 *}{origin}{origin of axes}{IN}
  2603. ** \param{Ppoint3 *}{xdirn}{x direction}{IN}
  2604. ** \param{Ppoint3 *}{y dirn}{y direction}{IN}
  2605. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  2606. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  2607. ** \paramend
  2608. ** \blurb{This function computes an object transformation which
  2609. ** maps unit vectors 
  2610. ** along the $x$, $y$ and $z$ axes onto unit vectors along the
  2611. ** corresponding axes 
  2612. ** of the new coordinate system.}
  2613. */
  2614. {
  2615.   Ppoint3 xaxis, yaxis, zaxis;
  2616.   Pmatrix3 objectxform;
  2617.  
  2618.   /* A new coordinate system is defined, centred at origin, its x-axis
  2619.   ** parallel to xdirn and the y-axis parallel to the component of ydirn
  2620.   ** which is perpendicular to the x-axis. The transformation matrix returned
  2621.   ** is an object transformation.
  2622.   */
  2623.  
  2624.   /* work out unit vectors along the axes of new coordinate system */
  2625.  
  2626.   xaxis = ptk_unitv3(xdirn);
  2627.   zaxis = ptk_crossv3(&xaxis, ydirn);
  2628.   zaxis = ptk_unitv3(&zaxis);
  2629.   yaxis = ptk_crossv3(&zaxis, &xaxis);
  2630.   yaxis = ptk_unitv3(&yaxis);
  2631.  
  2632.   /* form matrices - see "An Implementation of the GKS-3D/PHIGS Viewing
  2633.   ** Pipeline" for the maths. 
  2634.   */
  2635.  
  2636.   ptk_unitmatrix3(objectxform);
  2637.   objectxform[0][0] = xaxis.x;
  2638.   objectxform[1][0] = xaxis.y;
  2639.   objectxform[2][0] = xaxis.z;
  2640.   objectxform[0][1] = yaxis.x;
  2641.   objectxform[1][1] = yaxis.y;
  2642.   objectxform[2][1] = yaxis.z;
  2643.   objectxform[0][2] = zaxis.x;
  2644.   objectxform[1][2] = zaxis.y;
  2645.   objectxform[2][2] = zaxis.z;
  2646.   objectxform[0][3] = origin->x;
  2647.   objectxform[1][3] = origin->y;
  2648.   objectxform[2][3] = origin->z;
  2649.   ptk_concatenatematrix3(operation, objectxform, matrix, matrix);
  2650.  
  2651.   /* To output text, call above procedure to get the character plane 
  2652.   ** transformation:
  2653.   **
  2654.   **  ptk_0to3pt( text_point, text_direction_vector_1,
  2655.   **             text_direction_vector_2, PVReplacematrix,char_plane_xform)
  2656.   **
  2657.   ** As each point of the text is output, transform it by:
  2658.   **
  2659.   **  output_pipeline . char_plane_xform 
  2660.   */
  2661. }  /* ptk_0to3pt */
  2662.  
  2663. /*---------------------------------------------------------------------------*/
  2664.  
  2665. /*function:external*/
  2666. extern void ptk_oto3pt(C(Ppoint3 *) origin, C(Ppoint3 *) xdirn, 
  2667.              C(Ppoint3 *) ydirn, C(Pcompose_type) operation, C(Pmatrix3) matrix)
  2668. PreANSI(Ppoint3 *origin)
  2669. PreANSI(Ppoint3 *xdirn)
  2670. PreANSI(Ppoint3 *ydirn) 
  2671. PreANSI(Pcompose_type operation)
  2672. PreANSI(Pmatrix3 matrix)
  2673. /*
  2674. ** \parambegin
  2675. ** \param{Ppoint3 *}{origin}{origin of axes}{IN}
  2676. ** \param{Ppoint3 *}{xdirn}{x direction}{IN}
  2677. ** \param{Ppoint3 *}{y dirn}{y direction}{IN}
  2678. ** \param{Pcompose\_type}{operation}{concatenation operation}{IN}
  2679. ** \param{Pmatrix3}{matrix}{4x4 matrix}{OUT}
  2680. ** \paramend
  2681. ** \blurb{This function performs the same operation as
  2682. **  \pardesc{ptk\_0to3pt}, except the name has an \pardesc{o}\ 
  2683. ** (oh) instead of \pardesc{0}\ (zero). This function is provided for members
  2684. ** of the Fumbly Fingers Club.}
  2685. */
  2686. {
  2687.   ptk_0to3pt(origin, xdirn, ydirn, operation, matrix);
  2688. }  /* ptk_oto3pt */
  2689.  
  2690. /*--------------------------------------------------------------------------*/
  2691.  
  2692. /*function:external*/
  2693. extern void ptk_invertmatrix3(C(Pmatrix3) a, C(Pmatrix3) ainverse, 
  2694.                               C(Pint *)error)
  2695. PreANSI(Pmatrix3 a)
  2696. PreANSI(Pmatrix3 ainverse)
  2697. PreANSI(Pint *error)
  2698. /*
  2699. ** \parambegin
  2700. ** \param{Pmatrix3}{a}{4x4 matrix}{IN}
  2701. ** \param{Pmatrix3}{ainverse}{4x4 matrix}{OUT}
  2702. ** \param{Pint *}{error}{error code}{OUT}
  2703. ** \paramend
  2704. ** \blurb{This function computes the inverse of the $4 \times 4$
  2705. ** matrix \pardesc{a}, 
  2706. ** returning the result in \pardesc{ainverse}. 
  2707. ** If matrix \pardesc{a} is singular, then 
  2708. ** \pardesc{error} is set to $-1$, and \pardesc{ainverse} is undefined,
  2709. ** otherwise \pardesc{error} is set to \pardesc{ptkcpcok}.
  2710. ** This function uses Crout's method, a modification of Gaussian
  2711. ** elimination.}
  2712. */
  2713. {
  2714.   Pfloat b[4], c[4];
  2715.   Pfloat w, y;
  2716.   Pint z[4];
  2717.   Pint i, j, k, l, p;
  2718.  
  2719.   *error = -1;
  2720.   memcpy(ainverse, a, sizeof(Pmatrix3));
  2721.   for (j = 1; j <= 4; j++)
  2722.     z[j - 1] = j;
  2723.   for (i = 1; i <= 4; i++) 
  2724.   {
  2725.     k = i;
  2726.     y = ainverse[i - 1][i - 1];
  2727.     l = i - 1;
  2728.     p = i + 1;
  2729.     for (j = p; j <= 4; j++) 
  2730.     {
  2731.       w = ainverse[i - 1][j - 1];
  2732.       if (fabs(w) > fabs(y)) 
  2733.       {
  2734.     k = j;
  2735.     y = w;
  2736.       }
  2737.     }
  2738.     if (ptk_equal(y, 0.0))   /* matrix has no inverse */
  2739.       return;
  2740.     y = 1.0 / y;
  2741.     for (j = 0; j <= 3; j++) 
  2742.     {
  2743.       c[j] = ainverse[j][k - 1];
  2744.       ainverse[j][k - 1] = ainverse[j][i - 1];
  2745.       ainverse[j][i - 1] = -c[j] * y;
  2746.       ainverse[i - 1][j] *= y;
  2747.       b[j] = ainverse[i - 1][j];
  2748.     }
  2749.     ainverse[i - 1][i - 1] = y;
  2750.     j = z[i - 1];
  2751.     z[i - 1] = z[k - 1];
  2752.     z[k - 1] = j;
  2753.     for (k = 0; k < l; k++) 
  2754.     {
  2755.       for (j = 0; j < l; j++)
  2756.     ainverse[k][j] -= b[j] * c[k];
  2757.       for (j = p - 1; j <= 3; j++)
  2758.     ainverse[k][j] -= b[j] * c[k];
  2759.     }
  2760.     for (k = p - 1; k <= 3; k++) 
  2761.     {
  2762.       for (j = 0; j < l; j++)
  2763.     ainverse[k][j] -= b[j] * c[k];
  2764.       for (j = p - 1; j <= 3; j++)
  2765.     ainverse[k][j] -= b[j] * c[k];
  2766.     }
  2767.   }
  2768.   for (i = 0; i <= 3; i++) 
  2769.   {
  2770.     do 
  2771.     {
  2772.       k = z[i];
  2773.       if (k != i + 1) 
  2774.       {
  2775.     for (j = 0; j <= 3; j++) 
  2776.         {
  2777.       w = ainverse[i][j];
  2778.       ainverse[i][j] = ainverse[k - 1][j];
  2779.       ainverse[k - 1][j] = w;
  2780.     }
  2781.     p = z[i];
  2782.     z[i] = z[k - 1];
  2783.     z[k - 1] = p;
  2784.       }
  2785.     } while (k != i + 1);
  2786.   }
  2787.   *error = ptkcpcok;
  2788. }  /* ptk_invertmatrix3 */
  2789.  
  2790. /*--------------------------------------------------------------------------*/
  2791.  
  2792. /*function:external*/
  2793. extern void ptk_invertmatrix(C(Pmatrix) a, C(Pmatrix) ainverse, 
  2794.                               C(Pint *) error)
  2795. PreANSI(Pmatrix a)
  2796. PreANSI(Pmatrix ainverse)
  2797. PreANSI(Pint *error)
  2798. /*
  2799. ** \parambegin
  2800. ** \param{Pmatrix}{a}{3x3 matrix}{IN}
  2801. ** \param{Pmatrix}{ainverse}{3x3 matrix}{OUT}
  2802. ** \param{Pint *}{error}{error code}{OUT}
  2803. ** \paramend
  2804. ** \blurb{This function computes the inverse of the $3 \times 3$
  2805. ** matrix \pardesc{a}, 
  2806. ** returning the result in \pardesc{ainverse}. 
  2807. ** If matrix \pardesc{a} is singular, then 
  2808. ** \pardesc{error} is set to $-1$, and \pardesc{ainverse} is undefined,
  2809. ** otherwise \pardesc{error} is set to \pardesc{ptkcpcok}.
  2810. ** This function uses Crout's method, a modification of Gaussian
  2811. ** elimination.}
  2812. */
  2813. {
  2814.   Pfloat b[3], c[3];
  2815.   Pfloat w, y;
  2816.   Pint z[3];
  2817.   Pint i, j, k, l, p;
  2818.  
  2819.   *error = -1;
  2820.   memcpy(ainverse, a, sizeof(Pmatrix));
  2821.   for (j = 1; j <= 3; j++)
  2822.     z[j - 1] = j;
  2823.   for (i = 1; i <= 3; i++) 
  2824.   {
  2825.     k = i;
  2826.     y = ainverse[i - 1][i - 1];
  2827.     l = i - 1;
  2828.     p = i + 1;
  2829.     for (j = p; j <= 3; j++) 
  2830.     {
  2831.       w = ainverse[i - 1][j - 1];
  2832.       if (fabs(w) > fabs(y)) 
  2833.       {
  2834.     k = j;
  2835.     y = w;
  2836.       }
  2837.     }
  2838.     if (ptk_equal(y, 0.0))   /* matrix has no inverse */
  2839.       return;
  2840.     y = 1.0 / y;
  2841.     for (j = 0; j <= 2; j++) 
  2842.     {
  2843.       c[j] = ainverse[j][k - 1];
  2844.       ainverse[j][k - 1] = ainverse[j][i - 1];
  2845.       ainverse[j][i - 1] = -c[j] * y;
  2846.       ainverse[i - 1][j] *= y;
  2847.       b[j] = ainverse[i - 1][j];
  2848.     }
  2849.     ainverse[i - 1][i - 1] = y;
  2850.     j = z[i - 1];
  2851.     z[i - 1] = z[k - 1];
  2852.     z[k - 1] = j;
  2853.     for (k = 0; k < l; k++) 
  2854.     {
  2855.       for (j = 0; j < l; j++)
  2856.     ainverse[k][j] -= b[j] * c[k];
  2857.       for (j = p - 1; j <= 2; j++)
  2858.     ainverse[k][j] -= b[j] * c[k];
  2859.     }
  2860.     for (k = p - 1; k <= 2; k++) 
  2861.     {
  2862.       for (j = 0; j < l; j++)
  2863.     ainverse[k][j] -= b[j] * c[k];
  2864.       for (j = p - 1; j <= 2; j++)
  2865.     ainverse[k][j] -= b[j] * c[k];
  2866.     }
  2867.   }
  2868.   for (i = 0; i <= 2; i++) 
  2869.   {
  2870.     do 
  2871.     {
  2872.       k = z[i];
  2873.       if (k != i + 1) 
  2874.       {
  2875.     for (j = 0; j <= 2; j++) 
  2876.         {
  2877.       w = ainverse[i][j];
  2878.       ainverse[i][j] = ainverse[k - 1][j];
  2879.       ainverse[k - 1][j] = w;
  2880.     }
  2881.     p = z[i];
  2882.     z[i] = z[k - 1];
  2883.     z[k - 1] = p;
  2884.       }
  2885.     } while (k != i + 1);
  2886.   }
  2887.   *error = ptkcpcok;
  2888. }  /* ptk_invertmatrix */
  2889.  
  2890. /*--------------------------------------------------------------------------*/
  2891.  
  2892. /* end of tran.c */
  2893.